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ABSTRACT 

We  consider  the  nonlinear  programming  problem,  namely  minimizing  a 
nonlinear  function  subject  to  a  set  of  nonlinear  equality  and  inequality 
constraints.  Sequential  quadratic  programming  (SQP)  methods  are 
particularly  effective  for  solving  problems  of  this  nature.  It  is  assumed  that 
first  derivatives  of  the  objective  and  constraint  functions  are  available,  but 
that  second  derivatives  may  be  too  expensive  to  compute.  Instead,  the 
methods  typically  update  a  suitable  matrix  which  approximates  second 
derivative  information  at  each  iteration.  We  are  interested  in  developing  SQP 
methods  which  maintain  an  approximation  to  second  derivative  information 
projected  onto  the  tangent  space  of  the  constraints.  The  main  motivation  for 
our  work  is  that  only  the  projected  matrix  enters  into  the  optimality 
conditions  for  the  nonlinear  problem.  Updating  projected  second  derivative 
information  reduces  the  dimension  of  the  matrix  to  be  recurred;  we  avoid  the 
necessity  of  introducing  an  augmenting  term  which  can  lead  to  ill-conditioned 
matrices;  and  we  are  able  to  make  use  of  standard  quasi-Newton  updates 
which  maintain  hereditary  positive  definiteness.  We  discuss  four  possible 
formulations  of  the  quadratic  programming  subproblem  and  present 
numerical  results  which  indicate  that  our  methods  may  be  useful  in  practice. 
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Chapter  1 
Introductory  Material 


1.1  Introduction 

We  consider  the  problem  of  minimizing  a  nonlinear  function  subject  to  a  set  of 
nonlinear  constraints 

min  fix)  (NP) 

s.t.  Ci(x)  >  0       J  =  l,...,m; 

Cj(x)  =0       j  =  mi+\,  ■  ■  ■  ,m 
where  the  objective  and  constraint  functions  are  twice  continuously  differentiable. 

We  are  concerned  with  developing  algorithms  to  find  a  local  minimum  of  (NP)  by 

using  second  derivative  information  in  the  tangent  space  to  the  constraints. 

Sequential  quadratic  programming  methods  are  particularly  effective  for 
solving  problems  of  this  nature.  Algorithms  of  this  type  compute  the  minimizer  of 
(NP)  by  solving  a  sequence  of  quadratic  programming  subproblems.  It  is  assumed 
that  first  derivatives  of/  and  {c,}  are  available,  but  that  second  derivatives  may  be 
too  expensive  to  compute.  Instead,  the  methods  typically  update  a  suitable  n^n 
matrix  which  approximates  second  derivative  information  at  every  iteration. 

We  are  interested  m  developing  sequential  quadratic  programming  methods 
which  maintain  an  approximation  to  second  derivative  information  projected  onto 
the  tangent  space  of  the  constraints.  The  main  motivation  for  our  work  is  that  only 
the  projected  matrix  enters  into  the  optimality  conditions  for  the  nonlinear  problem. 
One    advantage    of    our    methods    is    that    updating    projected    second    derivative 


information  reduces  the  dimension  of  the  matrix  to  be  recurred.  Moreover,  we 
avoid  the  necessity  of  mtroducing  an  augmenting  term  which  can  lead  to  ill- 
.'onditioned  matrices  Since  the  reduced  matrix  of  second  derivatives  is  known  to  be 
positive  definite  at  the  solution,  we  are  able  to  make  use  of  standard  quasi-New  ton 
updates. 

Fou:  possible  formulations  of  the  quadratic  programming  subproblem  will  be 
presented  The  difficulties  arising  from  incomplete  knowledge  of  the  full  second 
derivative  matrix  will  be  discussed.  Proofs  of  global  convergence  of  some  of  these 
methods  for  a  restricted  class  of  problems  are  presented.  We  will  present  numerical 
result?  w'hi»_ii  indicate  that  our  algorithms  may  be  useful  in  practice. 


1.2  Notation 

Scalarh  will  usually  be  denoted  by  lower-case  Greek  letters,  e.g.    a,  p. 

Vectors  will  be  represented  by  lower-case  Roman  letters,  such  as  x,  y,  and  are 
assumed  to  be  column  vectors.  A  superscript  T  will  be  used  to  denote  transpose,  so 
x^  will  represent  a  row  vector.  Individual  components  of  a  vector  may  be 
referenced  by  the  subscripted  name  of  the  vector,  so  that  y^,  would  denote  the  third 
component  of  the  vector  y. 

I^/Iatrices  will  be  represented  by  capital  Roman  letters  A  particular  element 
wiil  be  referenced  by  the  name  of  the  matri.x,  m  either  upper  or  lower  case,  with 
tv.  ,ubscripts,  referring  to  the  row  and  column  in  which  the  element  appears,  e.g. 
<^i  i  '^r  ^i.i-   The  /th  row  of  the  matrix  A  will  be  denoted  by  a,. 

The  identity  matrix  of  order  n  will  be  denoted  by  /^  and  e,  will  be  used  to  refer 
to  its  /th  column;  R'^  will  be  used  to  represent  the  ^-dimensional  real  vector  space. 
Thf  scalai  €„^^  refers  to  the  machine  precision,  or  the  smallest  positive  value  on  a 
given  machine,  which,  when  added  to  one,  gives  a  result  greater  than  one. 
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We  will  make  use  of  the  following  vector  norms: 

X  111  =  Ski  .  (the  1 1  norm), 
1  =  1 

-^  II2  ~  (  ^^l  }^ '^  '  e^he  Euclidean,  or  2  norm), 
1  =  1 

X  Hoc  =    max  {  I  JT,-  I  }  ,  (the  °°  norm). 

ISi^n 

The  following  matrix  norms  will  be  used: 

^  II2  ~  (^raax('^'^'^))^^"  '  (the  largest  singular  value  of  A), 

-4  1!/^  =  (S  S<^n;)^'^  •  (the  Frobenius  norm). 
,  =  ij  =  i 

A  norm  referenced  without  a  subscript  will  be  assumed  to  mean  ||  .  II2. 

The  condition  number  of  a  matrix  A,  denoted  cond(A),  is  defined  as  ||y4||  ||A~^||. 

An    n^-n    matrix    A    is    positive  semi  — definite    if    x  Ax  "^  0    for    all    .vCR", 
positive  definite  if  x  Ax  >  0  for  all  nonzero  vSR". 

The  methods  that  we  will  be  studying  will  produce  a  sequence  of  points  in  R", 
X  ,x  ,  ,  converging  to  a  point  x' .    For  convenience,  we  will  assume  that  x  i'x* 

for  all  k.  Quantities  evaluated  at  iteration  /  will  be  denoted  by  a  superscript  of  i. 
The  following  definitons  of  convergence  refer  to  the  Q-rate  of  convergence  as 
defined  in  Ortega  and  Rheinboldt  (1970).  A  sequence  {.v*}  converging  to  x*  is  said 
to  converge 

linearly,  if 


W   ~  x^  II 

for  some  constant  0<ji,<l  and  all  k^k^; 
superlinearly ,  if 


II/  +  1  -  x*\\ 
lim  ".,  ,. j-f-  -  0; 


quadratic  ally ,  if 


k-'^     \\x     —  X 


^  M-, 


4- 


for  some  constant  ^.>0  and  all  k^kQ\ 
2 -step  superlinearly,  if 


fc-  k'~ -v 


-0: 


2-step  quadratically .  if 


■  \\- 
or  some  constant  fx>0  and  all  k'^k^. 

Given  any  function  h{x) .  we  use  the  notation  h  (x)^  or  Vh(x)  to  denote  the 
gradient  of  h  at  x. 

The  following  notation  will  be  used  in  reference  to  problem  (NP): 
Let  g{x)  =Vf(x),  the  gradient  of  the  objective  function, 
c(x)  =  [ci(i),  ■  ■  ,c„ix)]^.  A(.v)  =  [Vci(x)  •  ■  ,Vc„(x)]^,  the  matrix  whose  rows 
consist  of  ihe  transposed  constraint  gradients,  G(x)  =  V^/(.vi,  G,(.v)  =  V^c,(.r),  the 
Hessians  of  the  objecti\e  and  constraint  functions,  respectively.  We  will  assume  that 
the  unctions  /  and  c,  are  twice  differentiable  and  that  the  Hessians  of  the  objective 
an>.  constraint  functions  are  Lipschitz  continuous  matrix  functions  of  c. 

The    index    set   A//  =  {1,2,    ■■     ,mi}   refers   to   the   indices   of  the   inequality 
,:onstraints,  while  ME  =  {mi^\,  .  .m}  is   used  to  denote  the  equality  constraints. 


._et  A  be  partitioned  into 


-4/ 

A. 


,  where  A/  contains  the  transposed  gradients  of  the 


inequality  constraints  and  Aa-  denotes  the  transposed  gradients  of  the  equality 
constraints.  Similarly,  let  c;  and  c^  refer  to  the  values  of  the  inequality  and  equality 
constraints,  respectively. 

A  point  x  is  feasible  for  problem  <  NP)  if 

c/  >  0. 
Ce  =  0. 
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At  the  solution  to  (NP),  some  subset  of  the  inequahty  constraints,  as  well  as  all 
of  the  equality  constraints  are  zero.  These  constraints,  denoted  by  c(x  )  or  c  ,  are 
called  the  active  constraints  at  the  solution.  Any  algorithm  for  solving  (NP)  is 
concerned  with  identifying  the  active  set  at  the  solution.  Since  the  other  inequality 
constraints  have  a  positive  value,  they  remain  satisfied  within  some  neighborhood  of 
x' .  The  problem  then  essentially  reduces  to  the  problem  of  minimizing  the 
objective  function  with  respect  to  the  active  constraints.  The  index  set  J{x  ),  known 
as  the  working  set,  is  an  approximation  to  the  active  set  at  the  solution  which  is  used 
by  the  algorithm  at  the  point  a:*.  Any  constraint  which  is  in  the  working  set  is  called 
active,  although  it  may  not  actually  be  included  in  the  active  set  at  the  solution.  The 
constraints  in  the  working  set  are  referred  to  by  c(x  )  or  c  . 

Let  A{x)  denote  the  matrix  whose  rows  consist  of  the  transposed  gradients  of 
the  active  constraints  and  assume  that  A{x)  has  full  row  rank.  If  we  partition  Aj  into 


,   where  Aj  refers  to  the   inactive   constraint  gradients  and  Aj  denotes  the 


A  similar  partitioning  of  cj 


gradients  of  the  active  constraints,  then  A  = 

allows  us  to  express  c  as  {cj  c^)   . 

Let  Z{.v)  be  a  matrix  whose  columns  span  the  null  space  of  A(x),  and  let  Y{x) 
be  a  matrix  whose  columns  form  a  basis  for  the  range  space  of  A{x)^ .  Such 
matrices  can  be  obtained,  for  example,  by  forming  the  LQ  factorization  which 
produces 

A(x)Qix)  =  [L(x)  0].  (1.1) 

Let  r  be  the  number  of  active  constraints  and  let  t  =  n-r.    The  first  r  columns  of 

Q{x)  can  be  used  to  define  the  columns  of  Y{x),  while  the  columns  of  Z(x)  can  be 

formed  from  the  last  t  columns  of  Q{x). 

The  vector  gi(x)  and  the  matrix  W(j:,X)  will  be  used  to  denote  the  gradient  and 


-6- 

Hessian  of  the  Lagrangian  function,  which  will  be  introduced  in  the  next  section. 
The  approximation  to  the  Hessian  of  the  Lagrangian  function  at  r*  will  be  denoted 
by  H'' .  We  will  use  B''  to  refer  to  the  projected  Hessian  approximation. 


1.3  Optimality  Conditions 

The  following  presentation  of  the  optimality  conditions  refers  to  the  constraint 
qualifications  and  the  necessary  and  sufficient  conditions  for  a  point  x  to  be  a  local 
minimum  of  (NPi  as  presented  in  Fletcher  (1981)  We  will  assume  that  the  first 
and  second  order  constraint  qualifications  hold  at  x  .  A  sufficient,  but  not 
necessary,  condition  for  the  constraint  qualifications  to  hold  at  v  is  that  the 
gradients  of  the  active  constraints  are  linearly  independent,  i.e.,  A(x  )  is  of  full 
rank. 

A  point  -v*  is  optimal  for  iNP)  if 

I)      it  is  feasible 

2'      locally  there  is  no  other  feasible  point  with  a  lower  objective  function  value. 

Given  a  feasible  point  .v.  consider  a  step  along  a  first-order  feasible  direction  p. 
First,  consider  such  a  direction  which  maintains  to  first-order  a  zero  value  for  the 
constraints  in  the  working  set.   This  direction  must  then  satisfy 

A(x)p  =  0.  (1.2) 

The  matrix  Z{x)  forms  a  basis  for  all  p  that  satisfy  (1.2).  Such  a  direction  p  is  then 

of  the  form  p  =  Z(x)p^,  for  some  vector  p..    Consider  the  Taylor  series  expansion 

of/  at  x    along  such  a  direction  p: 

fix'  +  ap)  =  fix'  -h  aZ(.Y*)p,) 

=  fix')  +  ap]Zix*  f  gix' )  ^  Oia-p-)- 
In  order  for  x'  to  be  optimal, 


plZix')^gix*)  2:  0,    for  all  vectors  p.. 
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This  condition  reduces  to 


Zix'fgix')  =  0  (1.3) 

which  implies  that  g(x')  is  in  the  range  of  A{x')^ .  Therefore,  there  exists  a  vector 

X    such  that 

six')  =  A(x'f\' .  (1.4) 

The    quantity   Zix)^g{x)    is    called   the   projected  gradient   of  /   at   x.     Note   that 

requiring   the   projected   gradient   to   be   zero   is   a   special   case   of  the   optimality 

condition   ^(.v*)  =  0   for   unconstrained   optimization.     A    feasible   point   satisfying 

(1.3)  is  called  a  constrained  stationary  point. 

However,  it  may  be  that  the  current  working  set  is  not  the  active  set  at  the 
solution,  and  that  a  feasible  point  with  a  lower  function  value  could  be  found  by 
allowing  one  of  the  inequalities  in  the  current  working  set  to  take  on  a  positive 
value.    At  the  solution,  using  (1.4)  we  have 

XX,*a,(.v*)p  ^  0.  (1.5) 

A  first-order  feasible  direction/)  must  satisfy 

ai(x')p  =  0   idME, 
ai(x')p  >  0   z€A//. 
Condition  (1.5)  thus  implies 

X,   ^  0,     for  the  active  inequality  constraints. 
Therefore,  the  first-order  necessary  conditions  for  a  point  x     to  be  a  local 
minimum  of  NP  are  the  existence  of  a  vector  of  Lagrange  multipliers,  X  ,  such  that 

gix')  =  A(x'fK'  (1.6a) 

X,*  >  0,     for  the  active  inequality  constraints  (1.6b) 

c,(x')  >  0,  (1.6c) 

Ce(x*)  =  0.  (1.6d) 

A     pair     {x  ,  X  ),     satisfying     (1.6)     is     called     a    Kuhn  — Tucker  pair  of    (NP). 

Strict  complementarity  is   said  to  hold   if  there  is  no   ;   such  that  X,   =  c,(.r  )  =  0. 

Condition    (1.6a)    is   equivalent   to   requiring   x     to    be   a   stationary   point   of  the 


Lagrangian  function 

L(x,k)  ^  f(x)-X^iCi(x)  (1.7) 

when  \  =  X.  . 

Second-order  optimality  conditions  are  expressed  in  terms  of  the  Hessian  with 
respect  to  x  of  the  Lagrangian  function, 

W{xM  =  V^f(x)-2K^^c,{x). 
It  can  be  shown  that  a  necessary  condition  for  optimality  is  that  p^W(x' ,x')p  >  0 

for   any    feasible  p    (see    Fletcher    (198))).    Since  p  =  Z(x')p..   this    implies   that 

pJZix  )^W(x' ,k')Z{x*)p^  S:  0      for      all      p,,      i.e.,      the      projected      Hessian, 

Z(x'')^W(x' ,k')Z{x'),  is  positive  semi-definite.    Strict  positive  definiteness  of  the 

projected  Hessian  together  with  strict  inequality  in  (1.6b)  and  the  satisfaction  of  the 

first  order  conditions  is  a  second-order  sufficient  condition. 

In  addition  to  the  assumption  that  A  has  full  rank  (stated  earher),  we  shall 
assume  throughout  that  the  the  second-order  sufficient  conditions  hold  at  x  . 


1.4  Historical  Background 

The  earliest  approache';  to  solving  problem  (NP)  were  based  on  solving  a 
sequence  of  unconstrainc  •;  minimization  problems  whose  solution  in  the  limit 
approached  the  solution  of  (NP).  Penalty  function  and  barrier  function  methods 
belong  to  this  class  (see  Fiacco  and  McCormick  (1968j)  Penalty  function  methods 
minimize  a  sequence  of  subproblems  where  a  "penalty"  term  is  added  to  the 
objective  function  for  constraint  violation.  These  methods  generate  a  sequence  of 
infeasible  points.  The  penalty  term  grows  steadily,  forcing  convergence  to  a  feasible 
point  in  the  limit.  These  methods  are  not  appropriate  for  problems  in  which 
feasiblity  must  be  maintained  Barrier  function  methods,  on  the  other  hand, 
generate  a  sequence  of  feasible  points  by  creating  a  "barrier"  around  the  feasible 
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region.  Such  methods  involve  a  sequence  of  subproblems  where  the  function  to  be 
minimized  is  found  by  adding  a  weighted  sum  of  continuous  functions  with  a 
positive  singularity  at  the  boundary  of  the  feasible  region  to  the  objective  function. 
As  the  weight  assigned  to  the  barrier  term  is  decreased,  the  solutions  to  the 
unconstrained  subproblems  approach  the  constrained  minimum. 

The  availability  of  powerful  methods  for  unconstrained  optimization  and  the 
well-developed  theoretical  background  and  comparative  simplicity  of  these  methods 
is  attractive.  However,  in  practice  they  suffer  from  severe  numerical  difficulties, 
since  the  unconstrained  problems  become  increasingly  more  ill-conditioned  as  the 
solution  is  approached  (Murray  (1971),  Lootsma  (1972)). 

Algorithms  based  on  sequential  quadratic  programming  (SQP)  (also  known  as 
recursive  quadratic  programming  (RQP)  )  have  been  found  particularly  effective. 
Methods  of  this  type  approximate  the  minimizer  of  (NP)  by  solving  a  sequence  of 
quadratic  programming  (QP)  subproblems.  In  these  subproblems  a  local  model  of 
the  nonlinear  problem  is  posed  in  which  the  curvature  of  the  Lagrangian  function  is 
approximated  by  a  suitable  matrix  and  linear  approximations  are  made  to  the 
constraints.  The  solution  of  the  QP  is  used  to  provide  a  search  direction,  p*.  The 
algorithm  then  determines  a  step  length,  a*,  such  that  .v*  +  a*p*  is  a  sufficiently 
better  approximation  to  x  ,  in  a  sense  to  be  defined  later. 

At  the  solution  to  (NP)  we  are  concerned  only  with  the  constraints  that  are 
active.  Inactive  constraints  do  not  affect  the  solution  to  (NP),  since  they  remain 
satisfied  in  the  neighborhood  of  the  solution  and  thus  can  be  ignored.  Therefore, 
any  method  used  to  solve  (NP)  must  incorporate  some  means  of  determining  the 
active  set.  In  some  cases,  an  approximation  is  made  to  the  active  set  and  the  QP 
subproblem  is  posed  with  only  equality  constraints.  Such  subproblems  are  known  as 
EQPs  -  equality  constrained  quadratic  programs.  On  the  other  hand,  the  QP  can  be 
posed  as  an  IQP  -  an  inequality  constrained  quadratic  program.    Subproblems  of  this 
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nature  include  linear  approxraations  to  all  the  equality  and  inequality  constraints  of 
the  original  problem.  Although  the  terms  EQP  and  IQP  refer  to  the  formulation  of 
the  QP  subprobiem,  we  also  use  these  terms  in  a  looser  sense  to  mean  SQP  methods 
which  >o!ve  successive  QP  subproblems  of  either  type. 

The  EQP  .'■-ubproblem  for  (NP)  is  given  by 

mm    -p'^Hp  +  p'^g  (EQP) 

s.t.   A^p  =  -c. 

and  the  TQP  formulation  of  the  j:ubproblemi  is  given  by 

min    \p^Hp  +  p'^g  (IQP) 


s.t.  Alp  ^  -ci 
aIp  =  -ce 
where  H  is  some  approximation  to  the  Hessian  of  the  Lagrangian  function. 


Smce  the  Lagrangian  function  is  defined  in  terms  of  both  ,v  and  \,  these 
methods  must  incorporate  some  means  of  approximating  the  Lagrange  multipliers. 
At  the  solution,  the  Lagrange  multipliers  are  defined  by  (1.6a),  where  c(x*)  =  0. 
This  suggests  that  a  first-order  estimate,  X^.  can  be  found  at  .v*  by  solving  the  least 
squares  problem 

mxn\\A'\i-  g'Wl.  (1.8) 

The    least  squares  Lagrange  multiplier  estimates,  X.£,    are   only    an    estimate   of   the 
mulnplii-rs  defined  by  the  current  working  set. 

A  second  order  m  jltiplier  estimate,  •ti'^,  is  given  by  the  multipliers  defined  by 
the  solu'ion  p  of  the  QP  ~ubprob)em, 

A'\''  =  g'  +  Hp  (1.9) 

where  4  is  defined  by  the  working  set  at  the  solution  to  the  QP.    These  multiplier 

estimates  are  referred  to  as  the  QP  multipliers . 
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The  advantage  of  posing  the  subprobiem  as  an  EQP  is  that  the  solution  can  be 
easily  computed  by  solving  two  systems  of  linear  equations.  The  solution,  p  ,  can 
be  written  as 

p*  =  Ypy*  +  Zp* 
where  Y  and  Z  are  defined  as  in  (1.1).    The  vectors  py  and  p^  can  be  found  by 

solving 

Mp*  =  -c  (1.10a) 

Z^HZp'  =  -Z'^g  -  Z^HYp*  (1.10b) 

The  solutions  are  unique  when  A  has  full  row  rank  and  Z^HZ  is  positive  definite. 

This   method   of  computing  p*   is  discussed  by  Murray  and   Wright  (1982).    The 

solution  of  the  EQP  subprobiem  can  also  be  viewed  as  the  result  of  applying  one 

step  of  Newton's  method  to  the  nonlinear  system  of  equations 


g{x)-k{x)\ 

c{x) 
(see  Nocedal  and  Overton  (1985)). 


=  0  (1.11) 


In  general,  the  computation  of  the  solution  to  an  IQP  is  itself  a  finite  iterative 
process.  Essentially,  the  IQP  generates  a  sequence  of  working  sets  and  solves  the 
EQP  associated  with  each  working  set.  At  any  given  point,  a  search  direction  is 
generated  along  which  the  QP  objective  function  decreases  and  the  active  constraints 
remain  satisfied.  The  full  step  along  this  search  direction  may  cause  one  of  the 
currently  inactive  constraints  to  be  violated.  If  such  is  the  case,  a  step  is  taken  to 
the  closest  such  constraint.  The  newly-active  constraint  is  added  to  the  working  set 
and  the  process  is  repeated.  Eventually,  a  point  is  reached  which  is  the  minimum 
on  the  subspace  defined  by  the  current  working  set  (i.e.,  the  gradient  projected  onto 
the  tangent  space  to  the  constraints  in  the  working  set  is  zero).  To  determine 
whether  this  point  is  optimal  for  the  QP,  the  QP  Lagrange  multiplier  estimates  are 
checked.    If  some  multiplier  is  negative,  that  is  an  indication  that  the  QP  objective 
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function  can  be  further  decreased  by  deleting  the  corresponding  constraint  from  the 
working  set  This  process  of  revising  the  working  set  and  minimizing  along  the 
sub?  pace  defined  by  the  working  set  is  continued  until  a  point  is  reached  at  which 
the  jirojected  gradient  is  zero  and  the  Lagrange  multipliers  are  all  non-negative. 

£QP  methods  are  attractive  because  the  formulation  of  the  subproblem 
parallels  NewtonV  method  applied  to  an  unconstramed  minimization  problem. 
Moreover,  when  the  conect  working  set  is  identified,  Z(x' )'^W{x' .k')Zt,x')  must  be 
at  leaM  positive  semi-definite.  In  that  case  one  can  extend  quasi-Newton  updating 
methods  to  the  case  of  maintaining  a  positive  definite  approximation  to  Z^WZ. 
However,  the  difficulties  of  determining  the  correct  working  set  should  not  be 
minimi/.ed.  Algorithms  which  update  the  working  set  at  each  iteration  can  suffer 
from  repeated  insertion  and  deletion  of  the  same  constraint,  known  as  zigzagging 
(see  Wolfe  (1972)^ 

Methods  based  upon  solving  IQP  subproblems  are  not  dependent  on  a  strategy 
for  updating  the  working  set  at  each  iteration.  Such  algorithms  allow  for  fast 
recovery  from  a  bad  initial  approximation  to  the  active  set.  IQP  methods  typically 
involve  approximating  the  full  matrix  W,  since  the  constraints  defining  Z^WZ  change 
within  the  subprobleni.  However,  there  is  no  a  priori  reason  to  assume  that  \V  is 
positive  definite.  Moreover,  when  W  is  indefinite,  the  solution  may  not  be 
bounded,  and  the  direction  produced  by  the  QP  ma\  not  be  a  descent  direction. 
Therefore,  many  methods  which  update  the  tuU  second  derivative  matrix  of  the 
Lagrangian  function  add  a  term  of  the  form  pA^A  lo  the  definition  of  the 
Lagrangian  function  so  as  to  force  W  to  be  positive  definite.  This  augmenting  term 
can  lead  to  very  ill-conditioned  matrices,  although  in  practice  this  does  not  appear  to 
impede  the  performance  of  such  methods  (Giil,  ei.  al.,  private  communication). 

A  discussion  of  the  relative  merits  of  EQP  and  IQP  methods  can  be  found  in 
Vfv)rray   and   Wright   (1982)    and   Gil!   et.   al.   (1985cj.   In   the   neighborhood  of  a 
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solution,  one  would  expect  similar  performance  from  both  the  IQP  and  EQP 
methods.  Close  to  x'  it  is  expected  that  the  working  set  would  correspond  to  the 
correct  active  set.  Therefore,  the  solution  to  an  IQP  subproblem  would  be  found  by 
taking  one  step  along  the  subspace  defined  by  the  current  working  set  which  is 
equivalent  to  the  solution  of  the  EQP  defined  by  those  constraints. 

It  appears  that  SQP  methods  were  first  suggested  by  Wilson  (1963).  Wilson 
proposed  solving  an  IQP  at  each  iteration,  where  the  quadratic  approximation  to  the 
Lagrangian  function  was  formed  using  the  exact  Hessians  of  the  objective  and 
constraint  functions.  Murray  (1969),  Biggs  (1972)  and  Wright  (1976)  presented 
successive  EQP  methods.  Garcia  and  Mangasarian  (1976)  used  quasi-Newton 
updating  methods  to  approximate  the  full  Jacobian  of  (1.11).  However,  their 
approach  ignores  the  structure  of  the  Jacobian  matrix.  Moreover,  since  the 
Jacobian  is  always  indefinite,  the  DFP  and  BFGS  updates  cannot  be  used. 

Han  (1976,  1977)  proposed  a  successive  IQP  method  using  quasi-Newton 
updates  of  the  Hessian  of  the  Lagrangian.  However,  he  assumes  that  this  matrix  is 
everywhere  positive  definite.  Under  certain  conditions,  Han  proves  local 
superlinear  convergence  of  this  algorithm.  Powell  (1978)  presented  a  successive 
IQP  method  in  which  a  positive  definite  approximation  to  the  Hessian  of  the 
Lagrangian  is  maintained,  even  when  the  true  Hessian  is  indefinite.  This  method  can 
be  shown  to  have  local  two-step  superlinear  convergence  under  certain  conditions. 


1.5  Subproblems  Based  on  Projected  Hessian  Approximations 

We  are  interested  in  developing  algorithms  which  maintain  an  approximation 
to  the  projected  Hessian  rather  than  to  the  full  second  derivative  matrix.  Since  the 
projected  Hessian  must  be  at  least  positive  semi-definite  at  the  solution,  we  can 
make   use   of  standard  quasi-Newton   updates   which   maintain  hereditary  positive 
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dcfiniteness.  Updating  the  projected  Hessian  reduces  the  dimension  of  the  matrix  to 
be  recurred.  Moreover,  we  avoid  the  necessity  of  introducing  the  augmenting  term 
which  can  lead  to  ill-conditioned  matrices. 

The  idea  of  recurring  an  approximation  to  the  projected  Hessian  rather  than  to 
the  luil  second-derivative  matrix  was  suggested  by  Murray  and  Wright  f  1978).  Some 
partial 'ar  algorithms  and  their  convergence  analysis  were  given  for  the  equahty 
constrained  case  by  Coleman  and  Conn  (1984)  and  Noceda!  and  Overton  ( 1985). 

We  note  that  the  equality  constrained  case  is  much  simpler  than  the  general 
case,  since  the  active  set,  .vhich  determines  the  columns  of  Z,  is  known.  In  the 
inequalu}  constrained  case,  as  the  approximation  to  the  active  set  is  revised,  the 
dimension  of  Z^WZ  is  a!te?ed  Our  proposed  methods  are  designed  to  extend  the 
projected  Hessian  updating  algorithms  of  Nocedal  and  Overton  to  the  general 
inequality  constrained  problem 

In  'Uir  approach,  smce  ihe  Hessian  is  not  being  approximated  on  the  entire 
space  :ne  Z^HY  term  which  is  used  to  compute  p  in  EQP  methods  is  not  available. 
Similarly,  IQP  methods  require  computing  the  gradient  of  the  QP  objective  function 
at  intermediate  points,  which  is  obtained  by  multiplying  the  search  direction  by  the 
full  matrix  H ,  which  is  not  available. 

We  consider  severa!  ways  of  formulating  the  subproblem  to  accommodate  the 
fact  that  second  derivative  information  is  available  only  on  the  tangent  space  to  the 
constraints  Firstly,  we  can  formulate  the  subproblem  as  an  EQP.  replacing  Z^HZ 
by  the  approximating  term  S  and  ignoring  the  Z^HY  term.  The  drawback  of  this 
formulation  s  that  the  outer  iterations  must  incorporate  some  method  of  estimating 
the  active  set. 

On  the  other  extreme,  we  consider  posing  the  subproblem  as  an  IQP,  in  which 
the  full  Hessian  is  approximated  by  "unprojecting"  the  projected  Hessian,  i.e..  by 
defining  the  QP  as 
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1      ^ 

B    0 

min 

{p'Q 

piR- 

0   D. 

Q^p  +  p'^g 

s.t.   aJp  ^  -cj 

^Ip  =  -Ce 
where  B  is  an  approximation  to  the  projected  Hessian,  Q  is  defined  by  (1.1)  and  D 

is  a  positive  diagonal  matrix,  such  as  the  identity  matrix,  of  the  appropriate 
dimension.  Such  a  method  is  attractive  because  many  of  the  results  known  for  IQP 
methods  can  be  applied  to  the  special  case  where  the  positive  definite  matrix  is 
defined  as  above.  In  addition,  the  determination  of  the  working  set  is  done 
automatically  by  the  subproblem.  However,  since  the  matrix  is  not  accurate  on  the 
space  spanned  by  the  constraint  gradients,  the  search  direction  produced  by  the  IQP 
may  not  be  meaningful  when  far  from  the  solution.  Moreover,  the  selection  of  the 
working  set  is  predicated  upon  multiplier  estimates  which  may  not  be  reliable. 

A  major  problem  inherent  in  implementing  an  IQP  when  only  reduced 
Hessian  information  is  available  is  that  when  a  constraint  is  deleted,  it  is  not  clear 
how  to  expand  the  projected  Hessian.  When  the  full  Hessian  is  being  recurred,  an 
estimation  of  the  new  column  and  row  of  Z^HZ  can  be  obtained  by  expanding  the 
matrix  Z  and  multiplying  the  Hessian  onto  the  new  column  of  Z.  Since  we  do  not 
have  the  full  Hessian  available,  we  set  the  new  column  to  some  multiple  of  the 
corresponding  column  of  the  identity  matrix.  In  addition,  the  choice  of  constraint  to 
be  dropped  depends  upon  multiplier  estimates  which  may  not  be  reliable. 
Therefore,  we  consider  two  additional  approaches  which  are  midway  between  an 
EQP  and  an  IQP  method. 

In  the  first  method,  we  partially  solve  the  IQP  by  allowing  the  subproblem  to 
add  constraints,  but  terminate  the  iterative  process  as  soon  as  a  constraint  is  to  be 
dropped.  We  call  this  method  a  PQP  -  partial  QP  -  since  we  are  partially  solving  the 
IQP.  In  this  way  we  have  the  advantage  of  being  able  to  pick  up  many  constraints 
within   the  QP,  but  avoid   the  problems  of  expanding  the   projected  Hessian  and 
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relying  upon  multiplier  estimates  of  dubious  accuracy.    This  method,  however,  does 
entail  updating  the  active  set  approximation  outside  the  QP. 

In  the  second  method  we  allow  only  those  constraints  that  have  been  added 
during  the  current  QP  to  be  deleted.  Thus  the  approximating  matrix  is  allowed  to 
shrink  by  adding  constraints,  and  if  a  constraint  is  to  be  deleted,  second  derivative 
information  on  the  appropriate  subspace  can  he  recovered  from  the  projected 
Hessian  at  the  start  of  the  QP.  This  formulation  of  the  subproblem  will  be  called  an 
MQP  -  modified  QP.  Such  a  method  is  implemented  by  treating  all  constraints  that 
arc-  initially  active  as  equality  constraints,  while  including  the  inactive  constraints  as 
inequalities  Since  there  is  no  restriction  on  the  sign  of  a  multiplier  associated  with 
an  ^quality  constraint,  the  dropping  criterion  will  never  select  one  of  the  initially 
active  constraints  to  be  deleted 

In  the  neighborhood  of  a  solution,  we  would  expect  the  correct  active  set  to 
be  identified.  All  of  the  above  methods  would  then  find  the  solution  by  taking  one 
step  ii".  the  null  space  of  the  active  constraints.  Thus,  near  the  solution  all  four 
forniulations  reduce  to  an  EQP  method. 

Regardless  of  the  formulation  of  the  subproblem.  there  are  various  related 
issues  with  which  we  are  concerned.  In  later  sections  we  will  discuss  the  issues 
in'  olved  in  the  choice  of  update,  line  search  and  Lagrange  multiplier  estimates.  We 
w  ill  present  a  class  of  problems  which  provides  a  setting  in  which  the  performance 
of  our  algorithoiv  can  be  measured.  Convergence  properties  of  these  algorithms  will 
be  studied  Finally,  we  will  present  some  numerical  results  which  indicate  that  our 
algorithms  may  be  useful  in  practice. 
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Chapter  2 
Computational  Procedures 


2.1  Introduction 

NPSOL  is  an  SQP  software  package  developed  at  the  Systems  Optimization 
Laboratory  at  Stanford  University  (see  Gill  et.  al.,  1983a, b).  Given  an  initial 
estimate,  a:°,  to  the  solution  of  the  nonlinear  problem,  the  method  performs  a 
sequence  of  "major  iterations"  in  which  an  IQP  subproblem  is  solved  to  produce  a 
search  direction  p*.  The  new  estimate,  .r*^\  is  given  by  x*"*"^  =  x*  +  a*p*,  where 
a*  is  chosen  to  produce  "sufficient  decrease"  in  an  augmented  Lagrangian  merit 
function. 

The  QP  subproblem  itself  requires  an  iterative  procedure  to  solve  it,  and  thus 
several  "minor  iterations"  may  be  performed  within  the  QP  corresponding  to  a 
single  major  iteration.  The  matrix  H  used  in  the  QP  is  a  positive  definite 
approximation  to  the  Hessian  of  an  augmented  Lagrangian  function.  Evaluation  of 
the  objective  and  constraint  functions  and  updating  of  the  Hessian  are  performed 
only  in  the  outer  iterations. 

The  subroutine  package  QPSOL  used  to  solve  the  IQP  was  designed  for  use  in 
the  context  of  a  nonlinear  programming  algorithm,  and  therefore  includes  several 
features  tailored  to  the  problem.  Firstly,  the  linearization  of  the  nonlinear 
constraints  may  cause  the  QP  subproblem  to  be  infeasible,  although  there  may  be 
feasible  points  for  the  original  problem.  When  infeasibility  of  the  subproblem 
arises  from  the  nonlinear  constraints,  the  routine  temporarily  relaxes  the 
appropriate  constraints  in  order  to  allow  the  computation  to  proceed.    Secondly,  the 
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active  set  produced  b\  one  IQP  is  used  as  the  initial  active  set  for  the  next 
subproblem  In  the  neighborhood  of  a  solution  this  strategy  generally  allows  the 
IQP  to  be  solved  in  one  iteration.  In  addition,  the  QP  is  formulated  to  take  account 
of  the  special  structure  of  linear  and  bound  constraints. 

In  order  to  concentrate  on  the  issues  involved  in  restricting  the  second 
derivative  approximations  to  the  reduced  space  without  concerning  ourselves  with 
other  implementation  details,  we  have  developed  our  algorithms  by  modifying  the 
iMPS'.yi^  package.  We  have  endeavored  to  deviate  from  the  original  routines  as  little 
as  possible.  Since  our  algorithms  correspond  to  NPSOL  aside  from  the  differences 
induced  by  the  different  approximations  to  the  Hessian,  we  are  able  to  compare  the 
effect  of  'he  reduced  approximation  upon  the  performance  of  SQP  algorithms. 


2.2  The  TQ  Factorization 

A  crucial  part  of  many  quadratic  programming  methods  is  the  definition  of  a 
working  set  of  constraints  that  are  satisfied  exactly  at  the  current  iterate.  The 
workmg  set  typically  includes  a  mixture  of  bound  constraints  and  linear  constraints 
Changes  in  the  working  set  involve  changes  in  a  matrix  A  associated  with  the 
working  set.  In  general,  changes  in  the  matrix  A  are  reflected  by  updating  some 
factorization  of  .4. 

The   formulation   of  the  quadratic  programming  subproblem  used  in  NPSOL 
involves  constraints  of  the  form 

/  ^   ill)'  ^  "• 
The  matrix  A  is  therefore  composed  of  the  rows  of 

corresponding  to  the  active  constraints.    This  characterization  of  the  working  set 
enables  simple  hound  constraints  to  be  handled  separately.    However,  in  order  to 
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implement  such  methods,  the  factorization  must  be  maintained  subject  to  both  row 
and  column  updates.  The  TQ  factorization,  a  variant  of  the  standard  LQ 
factorization,  is  designed  to  enable  changes  in  the  working  set  to  be  implemented 
stably  and  efficiently. 

At  a  given  iteration,  the  working  set  consisting  of  r  constraints  will  include  a 
mixture  of  bound  and  general  constraints.  The  variables  corresponding  to  bound 
constraints  are  fixed  during  that  iteration,  while  the  remaining  variables  are  fi-ee. 
Let  the  subscripts  FX  and  FR  denote  the  fixed  and  free  variables,  respectively. 
Suppose  that  the  current  working  set  consists  of  rifx  bounds  and  m^  general 
constraints.  Let  nfjf  represent  the  number  of  free  variables,  and  let  M  denote  the 
matrix  whose  rows  are  the  active  general  constraints.  Since  the  variables  can  be 
reordered,  without  loss  of  generality  we  shall  assume  that  the  itfx  variables  are  the 
last  variables.  The  matrix  A  can  then  be  described  by 


A  = 


M 


M°    'S"      I-  (2.1) 


In  general,  SQP  methods  involve  the  use  of  a  matrix  Z  whose  columns  span  the 
null  space  of  A.  As  a  result  of  the  form  of  (2.1)  Z  may  be  expressed  in  terms  of 
only  those  columns  of  M  corresponding  to  the  free  variables.  Such  a  matrix  Z  is 
given  by 

2=\\'l  (2.2) 


0 
where  Zfj^  is  a  matrix  whose  columns  are  orthogonal  to  the  rows  of  Mfg. 

This  matrix  Z^^  is  obtained  from  the  TQ  factorization  of  Mf^,  which  is  defined 
by 

MfrQfr  =  (OT).  (2.3) 

The  matrix  2^/?   is  an  /i/r/jX^i^^  matrix,  and  T  is  an  mi'X-mi  reverse  triangular 

matrix,  i.e.,  7  is  a  lower-triangular  matrix  with  its  columns  in  reverse  order.    Note 
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chat  when  Q/r/f  is  orthogonal,  the  TQ  factorization  corresponds  to  a  permutation  of 
the  standard  LQ  factorization.  This  non-standard  approach  is  used  because  this 
Implementation  facilitates  updates  of  Z^^  when  a  general  constraint  is  deleted  from 
the  working  set  (see  Gill  et.  al.  1983c.  section  6.3). 

Let  nz  denote  n-r.  the  number  of  columns  of  Z.  It  follows  from  (2.2)  that  Z/r^ 
f  iefined  as  the  first  n^  columns  of  2/-/?  The  remaining  columns  of  Qf^  form  a 
basis  for  the  range  space  of  .A/Jr,  and  will  be  denoted  by  Yf/f 

fhe  TQ  factorization  i)f  '.he  full  matrix  .4  in  (2.1 )  is  given  by 


AQ 


0  Ifx  j        [Mfr  I^fx 


Z,,Y,,0]  0   0   /„ 

0    0  /„  0   T  M,J  ^       ' 


When  the  working  set  and  its  factorization  are  represented  as  above,  many  of 
the  calculations  involved  in  computing  the  search  direction  can  be  simplified.  Since 
the  searcf"  direction,  p.  is  defined  by 

Z^HZp,  =  -Z'^g 

P   =  Zp:, 

ii.  oollows  from  (2.2)  that  p  has  the  form  (pfK  0)'",  where  the  subscript  FR  is  used  to 
■  idicate  those  positions  of  the  original  vector  corresponding  to  the  free  variables. 
Also.  Z^g  =  Zfngfif  and  Z^HZ  =  ZjuHfuZfu.  Clearly,  the  amount  of  work 
involved  decreases  as  npx  increases.  Since  Z^HZ  is  generally  assumed  to  be  positive 
definite,  p^  is  computed  by  using  the  Cholesky  factorization 

Zi,H,f,Z,,  =  R^R.  (2.5) 

where  R  is  an  upper  triangular  matrix. 

As  constraints  are  added  to  and  deleted  from  the  working  set,  the  TQ  and 
Cholesky  factorizations  must  be  updated.  Adding  a  bound  constraint  to  the  working 
set  causes  one  column  of  Mf^  to  be  deleted  Therefore,  the  column  and  row 
dimensions  of  Zf^  as  well  as  the  dimensions  of  Qf^  and  R  are  decreased  by  one. 
The   dimension   of  T   is   unchanged,   although    the    matrix   is   updated.     Similarly, 
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deleting  a  bound  constraint  necessitates  adding  a  column  to  Mfr,  and  a  row  and  a 
column  to  Z^k  and  Qfr-  Adding  a  general  constraint  to  the  working  set  causes  the 
row  dimension  of  MpK,  the  dimension  of  T  and  the  column  dimension  of  Yfn  to 
increase  by  one  while  the  column  dimension  of  Z^^  decreases  by  one.  The  deletion 
of  a  general  constraint  reverses  these  changes. 

When  QpK  is  orthonormal,  sequences  of  plane  rotations  are  used  to  introduce 
zeros  into  the  appropriate  positions  of  these  matrices.  Any  rotations  applied  to  the 
columns  of  Yf^  are  also  applied  to  the  columns  of  7;  transformations  applied  to  the 
columns  of  Z/r/j  are  applied  to  the  columns  of  R.  The  matrix  Qf^  is  stored  in  a  large 
array,  ZY .  The  first  n^  columns  are  taken  to  be  Z^^,  and  K^/j  is  composed  of  the  last 
mi  columns.  Since  the  dimension  of  R  is  ^2  3"^  ^^e  dimension  of  T  is  mi,  both  R 
and  T  are  stored  in  a  single  array  RT .  The  first  nz  columns  are  used  to  store  R  and 
the  next  mi  columns  are  used  for  T .  The  advantage  of  this  implementation  is  that 
rotations  applied  to  the  columns  of  Q^n  can  be  applied  to  exactly  the  same  columns 
of  RT.  (Thus,  updates  to  the  columns  of  Zf^  will  affect  the  columns  of  R,  while 
updates  to  the  columns  of  Y/^^  will  be  applied  to  T .)  Detailed  descriptions  of  the 
updates  can  be  found  in  Gill  et.  al.  (1983c). 

An  option  included  in  NPSOL  involves  the  use  of 
stabilized  elementary  transformations  for  the  rotations  used  to  restore  T  to  triangular 
form.  (In  order  to  preserve  R^R,  the  updates  applied  to  restore  R  to  triangular 
form  must  be  performed  with  orthogonal  transformations.)  When  stabilized 
elementary  transformations  are  used,  Q  is  no  longer  orthogonal.  However,  in 
practice  it  generally  remains  fairly  well-conditioned. 


2.3  The  NPSOL  algorithm 

In  the  t\jllowing,  we  will  present  a  detailed  description  of  the  NPSOL  package 
>  Version  2.0,  1984). 

Th--;  problem  (NP)  is  treated  in  the  form: 

min  fix) 


subject  to   I 


Aix 
c(x) 


where /(x)  is  a  smooth  nonlinear  function,  A^  is  a  constant  matrix  representing  the 
iacobian  )f  the  linear  constraints,  and  c(x)  is  a  vector  of  smooth  nonlinear 
constraint    ^unctions.    Thus,    in    terms   of  the   notation   mtroduced    in   Chapter    1, 

..;  / 

c{x)  =  ■  A^x  ■  and  A(.r)  =      A^         The  user  must  supply  subroutines  that  define 

dx)  Aix) 

fix}.  c(x)  and  their  first  derivatives,  gix)  and  Aix),  as  well  as  x^,  an  initial  value 
for  x.  The  vectors  /  and  u  specify  lower  and  upper  bounds  for  the  variables  and 
-onstraints.  A  constraint  can  be  defined  as  an  equality  simply  by  setting  the 
corresponding  upper  bound  equal  to  the  lower  bound.  Elements  of  /  or  m  can  be  set 
to  special  values  that  will  be  treated  as  +^  or  —  oc  to  indicate  that  the  corresponding 
bounds  are  not  present.  .\  vector  of  feasibility  tolerances,  featol,  (supplied  by  the 
user)  is  used  to  determine  if  the  constraints  are  satisfied,  i.e..  if  the  ith  constraint  is 
violated  by  less  than /eafo/,,  it  will  be  considered  satisfied. 

We  use  a  superscript  to  denote  information  evaluated  at  each  major  iteration, 
and  a  subscript  for  quantities  modified  during  a  minor  iteration.  Thus,  for  example, 
A*  refers  to  the  matrix  of  active  constraint  gradients  at  the  (th  iteration  of  the  kth 
QP,  and  H  denotes  the  approximation  to  the  Hessian  matrix  at  the  kxh  major 
iteration,  which  is  not  changed  within  the  QP.  We  note  that  our  presentation  of 
NPSOL    is    not    intended    to    be    exhaustive.    For    a   thorough    description    of   the 
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implementation,  the  reader  is  referred  to  the  actual  Fortran  code. 
The  algorithm  proceeds  as  follows: 

1.  Starting  with  x^,  move  to  the  closest  point  that 

satisfies  the  linear  constraints.    This  move  defines 
the  initial  working  set,  J  . 

2.  Compute  the  TQ  factorization  of  A°,  where  the  rows 

of  A°  consist  of  the  gradients  of  the  constraints 
in  the  working  set,  J  . 

3.  Initialize  //°  by  setting  R°  =  I„. 

4.  Evaluate  f,  c^,  g^  and  A°. 

5.  Add  as  many  of  the  nonlinear  equality  constraints  and 

violated  nonlinear  constraints  as  possible  to  the 
working  set.   There  can  be  no  more  than  n  constraints 
in  the  working  set. 

FOR  /t  =  0  TO  maxiterations  DO 

6.  Compute  the  TQ  factorization  of  a'^ ,  obtaining 

y*.  Z*  and  r*.  Note  that  at  all  times,  A*, 
T\  Z*  and  K*  are  defined  by  7*. 
The  TQ  factorization  from  the  previous  iteration  can 
be  updated  by  refactoring  only  those  constraints  that 
have  changed  from  the  previous  iteration.   These 
constraints  include  all  the  nonlinear  constraints  in 
the  working  set,  as  well  as  linear  constraints  that  have 
changed  status. 


24 


7.         Solve  the  IQP  defined  by 


min  fnpip)  ^  TP^H'^p  +  g''\ 

piR'  ^ 

St.  I  -  c"  <  A^p  ^  u-c^ . 


where  c 


(See  Note  \ .) 


A^x*  \  and  A  *  = 


7.1  Compute  the  lower  and  upper  bounds,  /  and  /7, 

for  the  QP; 
Set    r=  l-c\ 
Set     .7  =  J/-C*. 

7.2  Step  10  the  constraints  in  the  working  set: 

Compute  r^  the  constraint  residuals,  defined  using 
the  appropriate  bound  for  each  constraint  in  /*. 

Solve  rV;=  -c^ 

Set    p  =  fVv'. 


(2.6) 


7  3         Find  the  nearest  feasible  point  for  the  QP  by 

solving  the  fallowing  linear  programming  problem: 

s.t.   Ij  ^  Ojp  ^  (~M  jiJ^ 
where  7i  is  the  set  of  indices  of  the  constraints 

violated  at  iheir  lower  bounds,  Jj  is  the  set  of 

indices  of  the  constrsints  violated  at  their  upper  bounds 

and  73  is  the  set  of  indices  of  the  satisfied 

constraints.  (See  Notes  1  and  2.) 

This  step  may  cause  additional  constraints  to  be 


(2.7) 
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added  to  the  working  set,  causing  7*,  and  thus 
A^  r*,  K*  and  Z*  to  be  modified. 

7.4  Set  /■  =  1. 

Define  j'l  =  ]\  Z\  =  Z^  Y^  =  K\  a\  =  A* 
and  T\  =  7*. 

7.5  Compute  the  Cholesky  factor.  R\,  of  Z'(h''Z\. 

(See  Note  3.) 

7.6  Set/g;^   =  fQp(p). 

Sqp,  =  SQpip)- 

where  the  objective  function  fgp  is 

defined  in  (2.6)  and  the  gradient  is  given  by 

gQpip)  ^  H^p  +  ^*.  (2.8) 

LOOP;       IF  Zfu'^l'l  is  positive  definite  and 
||Zf  gQp  II  is  sufficiently  small 

THEN 

7.7  Use  the  TQ  factorization  of  Af  to 
compute  the  QP  multipliers,  ti*,  defined  by 

A^Tl*  =  gQP,  (2.9) 

IF  they  are  all  sufficiently  large  and  have  the 
correct  sign  (see  Note  4) 
THEN  GOTO  EXITQP2. 
ELSE 

7.8  Delete  the  constraint  corresponding  to  the 
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multiplier  with  the  largest  violation  of 
its  sign  constraint  from  if. 
Update  the  TQ  factorization,  obtaining  an 
expanded  Zf,  and  expand  R^  accordingly. 
as  explained  in  Section  2.2. 
GOTO  LOOP 
ENDIF 


ENDIl 


7.9         IF  ||zf' i;j2/=J|  is  sufficiently  small  and 
zfH^-Z-  is  close  to  singular 
OR 
the  multipliers  •q'^  are  close  to  zero  and 
zfH'^Zi  is  close  to  singular 
OR 

i  >  maxiterations 
THEN  GOTO  EXITQPI. 
ENDIF 

■^.10        (Noi  used.) 

7.1 1  Compute  the  search  direction 

Solve  RfR-pi  =  -ZfgQp^. 
Set    pt  =  ZfPv 

7.12  Compute  the  maximum  step,  af,  along  this  direction. 

af  =  min(  i ,  9) 


where  9  =  min 


— T—         ajp,  <  0 
ajp, 

a^pf  >  0 
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7.13  Setp  =  p  +  afpf. 

7.14  StifQp^^fQpip) 

as  defined  in  (2.6)  and  (2.8). 

7.15  IF  af  =  9 

THEN 

Add  the  constraint  encountered  in  Step  7.12  to 

the  working  set,  giving  Vf+i- 
Update  the  TQ  factorization  to  get  Xf+i,  Zf+j, 
7"f+i  and  R-^^. 
(See  Note  5.) 
ENDIF 

7.16  Set  z  =  /+1 
GOTO  LOOP 

EXITQPl  -  set  error  code. 

7.17  EXITQP2  -  end  of  QP  subproblem. 

Define  p'^  =  p.    Set /^  =  J-. 

The  matrix  Af  is  defined  by  the  corresponding 

matrix  A*.   The  matrices  Tf,  Y^ 

and  Zf  are  defined  similarly. 

8.      Compute  the  value  of  the  slack  variables  that  will  minimize 
the  merit  function  Py  (see  Section  3.2). 
Compute  the  search  direction  for  the  slack  variables 
and  multipliers. 
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Couni  ihe  number  of  violated  constraints,  NVIOL, 
and  their  2-norm. 

9       Compute  |lZ*'^^*||,  !li?*||  and  the  least  squares  multipliers,  \*, 

defined  by  (1.8). 

10.  fF  IJZf'^^^il  ^  ice„^;,,  for  some  suitable  k,  and 

NVIOL  =  0  and 
the  mdhipliers  X*  .are  all  positive 
OR 
||zf^*|i  "small"  and  NVIOL  =  0  and 
the  multipliers  X^  are  all  positive 
and  the  previous  step  was  small 
THEN  GOTO  E.X!TNP2. 
BNDIF 

11.  Compute  the  maximum  step  d"  along  the  search  direction 

p^  subject  to  remaining  feasible  with  respect  to  the  bound 
and  linear  constraints, 

12.  (The  merit  function,  Pv(^''')  =  ^A/(a*,Ji•^^^P*).  will  be 

discussed  in  Section  3.2.) 
Compute  the  vectov  of  penalty  parameters,  p*.  such  that 


Compute  a*,  such  that  P^(jf*  +  oitV'')  is  sufficiently  lower 

than  P^(.v'')  and  a"  <  d'^. 
Set.v'=*^  =  x'  ^  a''p\ 
The  merit  function  evaluates /*'^',^*'^  ^c^^'^and  A*^". 
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If  an  acceptable  step  cannot  be  found,  GOTO  EXITNPl. 


13.     Update  the  approximation  to  the  Hessian  by  using  the 
BFGS  formula  (3.2),  where  y*  is  defined 
by  g*+l  -  gl,  the  change  in  the  gradient 
of  the  Lagrangian  function  using  the  QP  multipliers,  ^ 
If  necessary,  add  an  augmenting  term  to  keep  H 
positive  definite.   If  that  cannot  be  done,  the 
update  is  skipped. 


k 

k  +  i 


END  FOR 

EXITNPl  -  set  error  code. 
14.  EXITNP2  -  end  of  program. 

Notes: 

1.  In  fact,  constraint  i  is  considered  feasible  if  it  is  violated  by  no  more  than 

featoli- 

2.  If  no  feasible  pomt  can  be  found  for  (2.7).  that  could  be  due  to  inconsistent 
Hnear  constraints  in  the  QP.  In  that  case,  NPSOL  increases  featol  and  attempts 
to  solve  (2.7)  with  the  new  feasibility  tolerances.  This  process  is  repeated,  if 
necessary,  up  to  three  times.  If  a  feasible  point  was  still  not  found,  the 
program  exits  with  an  error  message. 

3.  The  projected  Hessian  is  formed  by  multiplying  out  zfH'z\,  where  Z\  is 
defined  by  the  current  working  set.  If  zfH'z\  is  not  positive  definite,  R\  is  set 
to  the  Cholesky  factor  of  the  largest  principal  positive  definite  submatrix  of 
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4.  The  optimality  condition  (1.6b)  restricts  the  multipliers  associated  with 
inequality  constraints  with  lower  bounds  to  be  nonnegative.  However,  the 
algorithm  considers  constraints  of  the  form  /,  ^  Cjix)  ^  w,,  which  essentially 
combmo  the  two  constraints  c,{x)  ^  /,  and  c,(x)  ^  «,.  Clearly,  the  upper 
bound  on  the  constraint  can  be  rewritten  as  -c,(.v)  >  =  -w,  Since  the 
multiplier  associated  with  this  reformulation  must  be  nonnegative,  it  follows 
that  the  muhiplier  associated  with  the  constraint  as  originally  written  as  an 
inequality  with  an  upper  bound  must  be  nonpositive.  Note  that  there  is  no 
restriction  on  the  signs  of  the  multipliers  associated  with  equality  constraints. 

5.  If  the  matrix  Z^  H'^Zf  defined  by  Zf  before  the  TQ  update  was  not  positive 
definite,  the  last  diagonal  element  or  the  whole  last  column  of  /?f  must  be 
recomputed. 


2.4  Description  of  the  Proposed  Algorithms 

In  the  following,  we  will  describe  our  four  proposed  algorithms  and  discuss 
how  they  are  implemented  in  the  context  of  NPSOL. 

The  IQP  Algorithm 

The  IQP  method  corresponds  most  closely  with  the  algorithm  used  in  NPSOL. 
The  major  difference  is  in  the  definition  of  the  positive  definite  matrix  used  in  the 
QP.  The  matrix  used  in  NPSOL  is  an  approximation  to  the  Hessian  of  the 
l.agrangian  function,  where  an  augmenting  term  is  added  if  the  Hessian  is  not 
positive  definite  Our  algorithm  maintains  an  approximation  to  the  Lagrangian  only 
in  the  tangent  space  to  the  constraint  gradients.    The  full  Hessian  is  approximated 
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by   expanding  the   projected  Hessian   into   an  nXn   matrix.     In   our   methods,   the 
update  is  applied  directly  to  the  projected  Hessian  approximation. 
Our  IQP  algorithm  is  described  below: 

1.  Initialize  W°  by  setting  /?°  and  Z°  =  /„. 

(See  Note  6.) 

2.  Follow  Step  1  of  NPSOL. 

3.  Follow  Step  2  of  NPSOL,  but  modify  /?°  as  well. 

(See  Note  6.) 

4.  Follow  Step  4  of  NPSOL. 

5.  Follow  Step  5  of  NPSOL,  but  modify  /?°  as  well. 

(See  Note  6.) 

FOR  )t  =  0  TO  maxiterations  DO 

6.  Follow  Step  6  of  NPSOL. 

7.  Solve  the  IQP  defined  by 


min/(2/>  =  ^P^H'^p  +  g^'p  (2.11) 

S.t.        I  -   c'^  <  a'^p   ^  u   -  c'^ 

where  A*  and  c*  are  defined  as  in  (2.6), 


H"  =  (2* 


B'     0 
0      D 


D  is  a  positive  diagonal  matrix  of  the  appropriate 
dimension,  Q'^  is  defined  in  (2.4)  and  fi*  is  the 
approximation  to  the  projected  Hessian  matrix. 
The  matrices  Q'^  and  S*  are  defined  by  the 


(2.12) 


-32- 

working  set  at  the  start  of  the  QP ,  7*.   The 
matrix  B'^  is  obtained  from  the  Cholesky  factors. 
R'^'r'^.    Within  the  QP,  changes  to  the  working  set 
will  affect  the  projected  Hessian  approximation 
and  the  TQ  factorization.  These  changes  are  made 
to  local  copies  of/?*,  2*  and  f*^,  denoted  hy 
P.  2*  and  f* 


7.1  -  7.2   Follow  Steps  7  1  -  7,2  of  NPSOL. 

7.3  Follow  Step  7.3  of  NPSOL.  but  modify  /?*  as  constraints  are 

added  to  the  working  set,  7*,,  (See  Note  6.) 

7.4  Follow  Step  ''A  of  NPSOL. 

1.5  Set  R\  =  R''. 

7.6-  7  16   Follow  Steps  7.6  -  7. 16  of  NPSOL. 

7.1''  Follow  Step   7.17  of  NPSOL.    The  Cholesky  factor  of  the 

projected  Hessian  approximation,  R  ,  is  set  to  /?, . 

8.  Follow  Step  9  of  NPSOL.    (See  Note  7.) 

9.  Follow  Step  8  of  NPSOL. 

10.  -  11.     Follow  .Steps  10  -  il  of  NPSOL. 

12.  Use  some  merit  function  lo  compute  a  step-length  a''. 

(The  choice  of  merit  function  will  be  discussed  in 
Section  3.2) 
If  an  acceptable  step  cannot  be  found,  GOTO  EXITNPl. 
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The  merit  function  evaluates  /^\  g''*'^,  c*""^  and  A*"'^ 

Set  7*^^  =  7^ 

13.  Compute  the  TQ  factorization  at  jr*^^  and 

use  some  updating  rule  to  obtain  /?*"^^  from  /?*. 
(The  updating  rules  will  be  discussed  in  Section  3.4) 

END  FOR. 

EXITNPI  -  set  error  code. 

14.  EXITNP2  -  end  of  program. 

Notes: 

6.  As  in  NPSOL,  /?*  is  used  to  provide  the  second-derivative  information. 
Hov^-ever,  in  the  NPSOL  algorithm  ^*  is  used  as  the  Cholesky  factor  of  the  full 
Hessian  approximation.  In  our  methods  /?*  contains  the  Cholesky  factor  of  the 
projected  Hessian  approximating  matrix.  Thus,  in  our  methods  ^*  must  be 
modified  as  described  in  Section  2.2  whenever  the  working  set  is  changed. 

7.  Steps  8  and  9  are  done  in  the  reverse  order  in  the  NPSOL  algorithm.  In  the 
NPSOL  algorithm,  the  search  direction  for  the  multipliers  is  obtained  from  the 
QP  multipliers,  ti*.  For  our  methods,  we  consider  using  the  least  squares 
multipliers  estimates  for  this  step.  Therefore,  we  first  compute  X.^  and  then 
compute  the  search  direction  for  the  multipliers. 

The  MQP  Algorithm 

Our  modified  QP  method  is  defined  similarly  to  our  IQP  method.  However, 
since  we  realize  that  our  approximation  to  the  full  Hessian  may  not  be  accurate  on 
the  full  space,  we  restrict  the  QP  by  not  allowing  it  to  drop  any  constraints  that 
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were  active  at  the  start  of  the  QP.  The  QP  may  add  constraints,  and  may  even 
delete  constramts  that  have  been  added  within  the  QP.  but  may  not  move  off  the 
constraints  that  were  active  previously  In  this  manner,  we  hope  to  minimize  the 
effect:'  of  dropping  constraints  based  upon  inaccurate  multiplier  estimates  and  of 
expanding  the  projected  Hessian  approximation  with  inaccurate  second  derivative 
information.  Constraining  the  QP  to  remain  on  the  original  active  set  is 
accumplished  by  treating  all  the  initially  active  constraints  as  equalities.  Since  there 
is  no  restriction  on  the  sign  of  the  multiplier  associated  with  an  equality  constraint, 
these  :oiistraint>  'a  ill  not  he  deleted 

We  note  that  it  :s  indeed  possible  for  the  QP  to  add  a  constraint  and 
5ubsequently  delete  it  This  scenario  occurs  when  there  is  another  constraint 
between  the  original  constraint  and  the  solution.  In  that  case,  the  second  constraint 
may  override  the  first.    'See  Figure  2.1). 

The  MQP  algorithm  is  described  as  follows: 

1.-5.    Follow  Steps  1  -  5  of  the  IQP  algorithm. 
FOR  k  =  0  TO  maxite rations  DO 

6.  Follow  Step  6  of  the  IQP  algorithm. 

7.  Solve  the  IQP  defined  by 

minfnp  (2.13) 

s.t.  Ojp  =--  I J  -  Cj.       jiL 

a]p  ^  Uj  -  c],     jdU 

where  /qp  is  defined  as  in  (2.1  l,i  and  W* 

is  defined  by  (2.12~i,    The  index  sets  L  and  U 

denote  the  constraints  that  are  in  the  working  set 
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at  the  start  of  the  QP,  Z*",  at  their  lower  and  upper 
bounds,  respectively.  The  index  set  /  refers  to  the 
inactive  constraints. 

7.1  -  7.17   Follow  Steps  7.1  -  7.17  of  the  IQP  algorithm. 

8.-11.      Follow  Steps  8  -  11  of  the  IQP  algorithm. 

12.  Use  some  merit  function  to  compute  a  step-length,  a  . 
The  merit  function  evaluates /■'\  g'''^\  c**^  and  A*""^. 

(See  Note  8.) 
Setx*"i  =  X*  +  aV*- 
Sety*"i  =  7^ 

13.  Compute  the  TQ  factorization,  projected  gradient  and 

constraint  values  at  the  new  point. 

14.  Compute  the  least-squares  multiplier  estimates.  If  the 

dropping  criteria  (see  Section  3.5)  are  satisified, 
drop  one  or  more  constraints  with  multipliers  of 
the  wrong  sign  (see  Note  4  of  the  NPSOL  algorithm) 
from  7*"^\  as  the  dropping  rule  specifies; 
update  the  TQ  factorization,  the  projected  gradient 
and  the  projected  Hessian  approximation. 

15.  Follow  Step  13  of  the  IQP  algorithm. 
END  FOR. 

EXITNPl  -  set  error  code. 

16.  EXITNP2  -  end  of  program. 
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Not«s: 

8.  Note  that  for  the  methods  which  do  not  solve  a  full  IQP  subproblem,  we  do  not 
stop  when  an  acceptable  steplength  cannot  be  found.  In  these  subproblems,  it 
nay  be  the  case,  for  example  that  the  QP  is  started  at  a  point  which  is  a  non- 
optima!  vertex,  and  terminates  immediately.  In  such  cases,  a  constraint  must 
le  deleted  from  the  working  set,  and  since  these  methods  do  not  delete 
-onstraints  within  the  QP,  the  working  set  must  be  modified  by  deleting  a 
ccnstraint  in  Step  14. 

The  PQP  Algorithm 

0- '  partial  QP  method  entails  posing  an  IQP  subproblem.  but  only  partially 
solving  '.t.  The  PQP  algorithm  avoids  the  problems  associated  with  dropping 
constraints  when  the  second  derivative  information  for  the  full  space  is  not 
ava.ilabli.  However,  by  allowing  the  PQP  to  partially  solve  the  full  IQP,  we  allow 
many  constraints  to  be  picked  up  in  the  course  of  one  major  iteration.  In  this 
manner  we  hope  to  retain  some  of  the  advantages  of  the  IQP  method,  while 
reduci'  g  the  possible  ill-effects  of  solving  the  full  IQP  when  the  full  Hessian  is  not 
being  approximated  The  PQP  is  implemented  by  following  the  steps  of  the  IQP 
until  .^dch  point  as  the  IQP  would  delete  a  constraint.  At  that  stage  the  PQP  is 
tex-minated. 

The  following  is  a  description  of  the  PQP  algorithm: 

1.-5.    Follow  Steps  1  -  5  of  the  IQP  algorithm. 
FOR  yt  =  0  TO  maxite rations  DO 

6.  Follow  Step  6  of  the  IQP  algorithm. 

7.  Partially  solve  the  IQP  defined  by  (2.11): 
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7.1  -  7,6      Follow  Steps  7.1  -  7.6  of  the  IQP  algorithm. 

The  partial  solving  of  the  IQP  is  accomplished  by: 

LOOP:  IF  llzf^e/'jl  '^  sufficiently  small 

THEN 

7.7  Compute  the  QP  multipliers,  t]'' ,  as  in  (2.9). 

7.8  GOTO  EXITQP2 
ENDIF 

7.9  -  7.15  Follow  Steps  7.9  -  7.15  of  the  IQP  algorithm. 

7.16  GOTO  LOOP   (See  Note  9.) 

EXITQPl  -  set  error  code. 

EXITQP2  -  Follow  Step  7.17  of  the  IQP  algorithm. 

8.-15.      Follow  Steps  8  -  15  of  the  MQP  algorithm. 

END  FOR 

EXITNPl  -  set  error  code. 
16.   EXITNP2  -  end  of  program. 

Notes: 

9.      If  the  full  step,  5*  =  1,  is  taken,  then  WzfligQpJ  =  0.  In  that  case  Steps  7.7 
and  7.8  will  be  executed  and  control  will  be  returned  to  the  outer  level. 

The  EQP  Algorithm 

For  our  EQP  method,  we  simply  solve  the  EQP  defined  by  the  current  working 
set  and  ignore  the  Z^WY  term  (as  described  in  Section   1.4).  This  algorithm  then 
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amounts  to  taking  one  step  of  an  IQP.  The  justification  for  ignoring  the  Z^'WY  term 
is  that  near  the  solution  we  expect  the  active  constraints  to  be  satisfied,  causing  the 
Py  ten  •  to  be  near  zero. 

A   'escription  of  the  fc"QP  algorithm  is  presented: 

1  -  5.     Follow  Steps  1  -  5  of  the  IQP  algorithm. 
FOR  ;^  =  0  TO  inaxhf rations  DO 

6.  Follow  Step  6  of  the  IQP  algorithm, 

7.  Solve  the  F.QP  which  is  defined  by  performing  one  iteration 

of  the  TQP  defined  by  (2.11): 

7.1  -  7.3      Follow  Steps  7,1  -  7  3  of  the  IQP  algorithm, 

7.4  Follow  Step  '!  6  of  the  IQP  algorithm,    (See  Note  10.) 

7  5  Solve 

R''r'p^=  -Z'^gQp  (2.14) 

Set  ^*  =  Z'^p^ 

7.6  Compute  «*■.  the  maximum  step  along  p*^. 

as  defined  in  Step  '',  12  of  the  IQP  algorithm. 

7.7  Setp*  =  p'^  +  aV'' 

7.8  Follow  Step  7, 14  of  the  IQP  algorithm , 

7.9  IF  a*  #  1 

THEN 

Add  the  new  constraint  to  the  working  set,  7". 
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Update  7*,  (2*  and  /?^ 
ENDIF 

7.10  Compute  ti*^,  the  QP  multipliers  as  defined  in  (2.9). 

7.11  EXITQP2  -  Follow  Step  7.17  of  the  IQP  algorithm. 
8.  -  15.       Follow  Steps  8  -  15  of  the  MQP  algorithm. 

END  FOR. 
16.     EXITNP2  -  end  of  program. 

Notes: 

10.    The  EQP  subproblem  is  solved  in  one  step.    Therefore,  there  is  no  need  for 
subscripts  on  the  QP  quantities  to  indicate  the  minor  iteration  number. 


2.5  Example 

In  order  to  compare  and  contrast  these  four  methods,  let  us  examine  the 
behavior  of  each  of  these  algorithms  on  a  small  example.  Consider  the  following 
problem: 

min  f(x)=xj  -i-  xj  +  x] 


xiH^ 

s.t. 

l<Xi 

1<a:2 

<4 

^3 

<4 

0<.<:i 

-  ^2 

0<xi 

v2 
^2 

with  .v°  =  (4,  2,  4)^".   The  solution,  x*,  is  at  (1,  1,  0)^. 
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rhe  IQP  Method 

At  the  start  of  the  algorithm  (Step  2  of  the  IQP)  all  the  linear  and  bound 
constraints  which  are  initially'  on  their  bounds  are  added  to  the  working  set. 
Therefore,  the  bound  constrsint  x;,  =  4  is  picked  up  immediately. 

The  first  QP  subprobiem  is  posed  with  .13  =  4  initially  active.  Since  any  non- 
zero step  along  the  negative  projected  gradient  ((  — 3;cr,  -Z.vi,  0)  )  would  violate 
the  nonlinear  constraint  Vj  ^  xl,  that  linearized  constraint  is  immediately  added  to 
ithe  working  set.  Keeping  the  Hneanzed  constraints  active,  the  QP  steps  along  the 
negative  gradient  in  the  tangent  space  to  the  constraints  and  encounters  the 
linearized  constraint  for  vj  =  X2  at  the  point  (1.333,  1.333.  4)^".  With  this  constraint 
added  to  ihc  working  set,  the  number  of  active  constraints  equals  the  dimension  of 
the  problem  and  the  projected  gradient  is  zero.  However,  the  QP  multipliers  are 
dot  all  positive  and  the  constraint  -.]  =  xf  is  dropped.  The  QP  now  moves  along  the 
lines  .'.  =  .V;  and  X3  =  4  until  the  constraint  x-.  —  I  is  reached.  Once  again,  there 
are  three  active  constraints.  .A.t  this  point  the  bound  .V3  =  4  is  deleted  from  the 
working  set  and  the  QP  moves  to  ( 1 ,  1,  -4)^.  The  projected  gradient  at  that  point 
is  .'.ero  and  all  the  QP  multipliers  are  positive,  so  the  QP  is  terminated. 

The  full  Newton  step  to  r  =  (I,  1,  -4)^  is  accepted  by  the  line  search  in  Step 
12.  Although  the  constraints  are  all  satisfied,  the  projected  gradient  is  still  "large". 
The  second  QP  is  posed  with  the  constraints  X]  =  X2  and  v;  =  1  active,  and  steps  to 
(1,  1,0).  The  outer  iteration  again  accepts  the  full  step  to  jr^  =  t  1,  1.  0)^,  which  is 
the  solution.  In  general,  the  method  does  not  produce  the  exact  solution,  but 
generates  an  infintte  sequence  oi  iterates  which  converge  to  x  . 

The  EQP  Method 

While  the  iQP  method  may  make  many  changes  to  the  working  set  at  each 
iteration,  the  EQP  method  can  take  only  step  for  each  QP.    That  step  may  be  a  full 
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Newton  step  in  the  tangent  space  to  the  active  constraints,  or  it  may  be  a  step  to 
another  constraint.  In  the  EQP  method,  constraints  are  dropped  outside  the  QP. 

As  in  the  IQP,  the  bound  X3  =  4  is  piciced  up  at  Step  2.  The  first  QP 
subproblem  immediately  encounters  the  nonlinear  constraint  xi  =  xl-  However, 
unlike  the  IQP,  the  QP  terminates  after  this  one  step,  giving  x^  =  x^  =  (4,  2,  4)^. 

Although  the  multipliers  are  not  all  positive,  the  projected  gradient  is  quite 
large,  causing  the  dropping  criteria  of  Step  14  not  to  be  satisfied.  The  constraints 
Xl  =  X2  and  .V3  =  4  are  used  to  define  the  search  direction  in  the  QP.  The  QP  steps 
along  the  negative  projected  gradient  in  the  tangent  space  to  those  constraints,  and 
encounters  the  constraint  x^  =  jtt  at  the  point  (1.333,  1.333,  4)^.  Again,  the  EQP 
ends  after  this  one  step,  and  x^  =  (1.333,  1.333,  4)^. 

Since  there  are  now  three  constraints  in  the  working  set  and  the  dimension  of 
the  problem  is  three,  the  projected  gradient  is  zero.  Therefore,  before  posing  the 
next  QP  subproblem,  the  EQP  method  checks  the  least  squares  multiphers  in  order 
to  determine  which  constraint  should  be  dropped  from  the  working  set  and 
accordingly  drops  the  bound  constraint.  The  third  QP  is  thus  posed  subject  to  the 
two  constraints  in  the  working  set,  JCj  =  Xj  and  Xj  =  xi-  Since  the  constraints  in  the 
working  set  are  not  exactly  satisfied  (.tj  ¥=  xl),  a  step  is  taken  to  the  linearized 
constraints  (Step  7.1).  The  second  component  of  the  QP  step  is  taken  to  minimize 
the  quadratic  approximation  to  the  objective  function  in  the  tangent  space  to  the 
active  constraints,  which,  in  this  case,  is  essentially  a  step  along  x^.  The  full  QP 
step  is  taken  to  give  .v^  =  (1.0667,  1.0667,  -4)^. 

Once  again,  the  projected  gradient  is  too  large  to  allow  for  deletions  from  the 
active  set,  and  the  next  QP  performs  similarly.  A  step  is  taken  to  the  linearized 
constraints,  followed  by  a  step  along  .V3,  providing  the  new  iterate 
x"^  =  (1.0039,  1.0039,  0)^. 
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Tbe  following  two  subproblems  step  to  the  linearized  constraints,  producing  the 
iterates  x^  =  (1.00015,  1.00015,  0)^  and  x^  -  (1.0000.  1.0000,  0)^.  At  this  point,  all 
the  constraints  are  satisfied  and  the  projected  gradient  is  zero.  The  least  squares 
multipliei  estimates  indicate  that  the  nonlinear  constraint  cj  =  xl  can  be  deleted 
from  the  working  set.  The  next  QP  attempts  to  take  a  step  in  the  new  space  defined 
by  the  new  working  set,  but  immediately  encounters  the  bound  constraint  Vj  =  1. 
Picking  up  the  bound  causes  the  projected  gradient  to  be  zero  again,  and  since  all 
the  least  squares  multiplier  estimates  are  now  positive,  the  program  terminates  with 
.r'  =  (1,  1,0)^. 

Thf  POP  Method 

As  m  the  other  cases,  the  PQP  starts  with  the  bound  x^  =  4  initially  active. 
The  first  QP  immediately  adds  the  nonlinear  constraint  xj  =  xj-  However,  unlike 
the  EQP.  the  PQP  does  not  terminate  after  this  one  action,  but  continues  to  move 
along  the  linearized  constraints  so  as  to  minimize  the  quadratic  objective  function. 
This  path  is  followed  until  it  reaches  the  constraint  Xj  =  X2.  which  is  then  added  to 
the  working  set.  Since  three  constraints  are  now  active,  the  projected  gradient  is 
zero  However,  unlike  the  IQP,  the  PQP  may  not  delete  any  constraints  within  the 
subproblem,  and  therefore  terminates.  The  new  iterate  is  thus  given  by 
x^  =  (1..333,  1..333.  4)" 

Before  posing  the  second  subproblem,  the  PQP  checks  the  multiplier  estimates 
and  drops  the  bound  x^  -  4.  .A.s  did  the  EQP,  the  PQP  first  steps  to  the  constraints 
and  then  along  the  negative  gradient  in  the  tangent  space  to  the  active  set,  to 
provide  .<:^  =  (1.0667.  1.0667,  -4)^. 

Since  no  new  constraints  are  picked  up  until  the  last  step,  the  PQP  proceeds 
identically  to  the  EQP,  producing 

x^  =  (1.0039,  ).0O39,0)^,  r"*  =  (1.00015,  1.00015,0)^  and.  finally, 
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x^  =  v*  =  (1,  1,0)^. 

The  MQP  Method 

At  each  iteration,  the  MQP  method  solves  a  subproblem  in  which  the 
constraints  in  the  working  set  are  included  as  equalities  and  the  remaining 
constraints  of  the  original  problem  are  included  as  inequalities. 

The  bound  X3  =  4  is  initially  active,  as  was  the  case  with  the  other  methods. 
The  first  QP  proceeds  as  in  the  IQP  case,  by  picking  up  the  nonlinear  constraint 
Xi  =  X2  and  moving  along  the  linearization  of  that  constraint  until  it  encounters  the 
linear  constraint  Xi  =  xj  at  the  point  (1.333,1.333,4)'^.  Since  the  projected 
gradient  is  zero  and  all  the  QP  multiplier  estimates  are  not  positive,  the  MQP  drops 
the  constraint  x^  =  xl-  This  constraint  may  be  dropped  since  it  was  not  in  the 
working  set  at  the  start  of  the  QP.  As  did  the  IQP,  the  MQP  steps  along  ;ci  =  xj 
and  x^  =  4  until  the  bound  Xj  =  1  is  reached.  Once  again,  the  projected  gradient  is 
zero  and  there  is  a  negative  QP  multiplier  estimate.  However,  although  the  IQP 
would  now  delete  x^  =  4,  the  MQP  cannot,  since  that  bound  was  active  at  the  start 
of  the  QP.   The  MQP  thus  terminates  and  the  new  iterate,. v^  is  given  by  (1,  1,  4)^. 

Before  posing  the  second  subproblem,  the  method  checks  the  projected 
gradient  and  the  multiplier  estimates,  to  determine  whether  constraints  should  be 
deleted  from  the  working  set.  At  this  point,  the  bound  X3  =  4  is  deleted.  The 
second  QP  moves  along  .1:3  to  the  point  (1,1,  -4)^.  However,  although  the  full  step 
defined  by  the  QP  had  been  accepted  by  the  line  search  in  all  the  other  cases,  at  this 
point  the  constraint  and  objective  function  values  are  identical  at  the  points 
(1,  1,4)  and  (1,  1,  — 4)  .  Thus,  the  full  step  does  not  provide  a  "sufficient 
decrease"  in  the  merit  function.  The  line  search  backtracks  along  the  QP  search 
direction,  and  accepts  the  point  halfway  along  the  full  step,  moving  to 
x    =(1,1,  0)^,  which  is  the  solution. 
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Summai  y 

In  examining  the  behavior  of  these  four  methods  on  the  above  example,  we  are 
able  to  -ee  that  the  various  methods  actually  do  perform  differently  in  practice.  In 
particubii  the  example  demonstrates  how  an  incorrect  approximation  to  the 
working  'Ct  can  cause  the  EQP  to  converge  >>!ow!y,  and  the  advantage  of  the  IQP  in 
being  able  )o  quickly  determine  the  correct  active  set.  However,  this  example  is  not 
representative  of  the  general  case  In  more  highly  nonlinear  problems,  the  IQP 
method  may  suffer  from  revising  the  active  set  based  on  the  linearizations  to  the 
constraints  which  may  be  very  inaccurate  when  far  from  the  solution  A  discussion 
of  the  advantages  and  disadvantages  of  both  IQP  and  EQP  methods  can  be  found  in 
Murray  and  Wright  (1982)  and  Gill  et.  al.  (1985c). 
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Figure  2.1.  This  picture  depicts  the  difference  between  the  directions  produced  by 
the  MQP  and  PQP  methods  The  contours  of  f^p  are  shown,  along  with  two  linear 
constraints.  The  unconstrained  solution  would  be  at  a.  The  solution  to  the  PQP  is 
found  at  the  intersection  of  the  two  constraints,  b.  For  the  MQP  method,  the  algo- 
rithm can  move  off  the  first  constraint  and  find  a  better  solution  at  c. 
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Chapter  3 
Implementation  Issues 


3.1  Multiplier  Estimates 

Lagrange  multiplier  estimates  are  used  to  update  the  Hessian  approximation, 
to  determine  whether  constraints  should  be  deleted  from  the  active  set  and  as 
parameters  of  the  augmented  Lagrangian  penalty  function.  Since  the  definition  of 
the  objective  function  of  the  subproblem  depends  upon  the  second  derivative  matrix, 
convergence  of  SQP  methods  depends  critically  upon  convergence  of  the  multiplier 
estimates. 

At  X  ,  the  Lagrange  multipliers  are  defined  by 

Aix'fi'  =  g{x') 
where 

cix*)  =  0  . 
This  suggests  that  a  first-order  estimate,  X^,  at  .v*  can  be  found  by  solving  the  least 

squares  problem  (1.8).    This  system  is  compatible  at  any  point  where  the  gradient  is 

a  linear  combination  of  the  rows  of  a''.   Therefore,  even  when  the  residual  is  zero, 

X.*  is  only  an  estimate  of  the  multipliers  defined  by  the  current  active  set,  unless 

cix'')  =  0.    The  TQ  factorization  of  A  can  be  used  to  compute  X*.    If  A(.v*)  has  full 

rank  and  ||a:*  -  x'  ||  is  of  order  e  for  sufficiently  small  €,  then  ||  X^  -  X*  ||  is  of 

order  €.    This  estimate  is,  thus,  a.  first-order  estimate  of  X  . 

A  second- order  multiplier  estimate,  t]'' ,  is  given  by  the  multipliers  defined  by 
the  solution  p  of  the  QP  subproblem,  as  in  ( 1 .9) .    If  A*  is  of  full  rank,  .t*  and  X*  are 
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sufticiently  close  to  x    and  \  ,  respectively,  and  p'^  is  a  second-order  estimate  of  the 
step  from  v'^  to  v  .  t|*  will  be  a  second-order  estimate  of  X*  (see  Gill  et.  al.  (1981)). 

A  second-order  multiplier  estimate  at  a*  is  a  prediction  of  the  first-order 
multipliers  at  the  new  point.  It  should  be  emphasized  that  t]*'  depends  critically 
upon  p  and  for  a  poor  choice  of  p''  can  be  complete!)  maccurate.  Moreover,  when 
the  projected  Hessian  is  recurred,  second  derivative  information  is  not  available  for 
the  enDre  space  and  thus  the  second-oraer  multiplier  estimate  may  lot  be  reliable 
Wrigh;  (1976)  shovvs  that  using  the  least  squares  multiplier  estimates  at  the  new 
point  for  updating  the  Hessi.ir.  approxunation  cannot  impede  the  faster  rate  of 
convergence  of  the  underlying  method  Using  a  different  approach,  Goodman 
( 1985'  also  shows  that  a  \'ewton-type  method  thai  uses  X^  is  quadratically 
convergent.   Therefore,  we  find  it  preferable  to  use  first-order  multipliers  at  ,v*"^^ 


i.2  Merit  Functions 

An  essential  element  of  any  algorithm  for  solving  optimization  problems  is  a 
technique  for  forcing  convergence  from  a  poor  starting  point.  Thi>-  is  generally 
accomplished  3y  viewing  the  step  to  the  new  point  as  a  search  direction  and 
choosing  a  step  length  along  that  path  in  such  a  manner  so  that  the  new  point  is  a 
"sufficiently  better"  approximation  to  the  solution.  In  the  unconstrained  case, 
progress  towards  a  solution  can  be  measured  in  terms  of  the  decrease  in  the 
objective  function.  The  new  iterate  can  therefore  be  chosen  as  any  point  along  the 
line  where  an  "acceptable"  reduction  in  the  obiective  function  is  attained.  Linearly 
constrained  problems  are  generally  solved  by  computing  an  initial  feasible  point  and 
geneiating  subsequent  iterates  so  as  to  continue  to  satisfy  the  con.straints  Thus, 
only  the  objective  function  needs  to  be  considered  in  computing  the  step  length. 
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When  nonlinear  constraints  are  present,  however,  it  is  not  straightforward  to 
produce  a  sequence  of  points  that  remain  feasible.  The  line  search  routine  must 
therefore  balance  the  reduction  of  the  objective  function  against  the  satisfaction  of 
the  constraints.  Various  merit  functions  have  been  proposed  for  measuring  the 
quality  of  a  step  length.  Ideally,  the  merit  function  should  be  chosen  so  as  to 
enforce  global  convergence.  It  is  desirable  that  it  be  easy  to  implement  and  not  too 
costly  to  evaluate.  At  each  iteration,  the  direction  produced  by  the  subproblem 
must  be  a  descent  direction  for  the  merit  function.  The  behavior  of  the  merit 
function  should  not  inhibit  the  convergence  of  the  underlying  method.  Asymptotic 
convergence  rates  are  based  upon  the  assumption  that  in  the  neighborhood  of  a 
solution  step  lengths  of  unity  are  taken.  It  is  to  be  hoped  that  near  a  solution  taking 
a  step  of  length  one  will,  in  fact,  result  in  a  sufficient  decrease  in  the  merit  function. 

The  quadratic  penalty  function  (Courant  (1943))  is  defined  by 

PQix,p)  ^f(x)  +  pc,ixfc,ix) 
where     c^.(x)     are     the     violated     constraints.      This     function     is     continuously 

differentiable,  but  has  discontinuities  in  the  second  derivative  at  points  where  a 

constraint  becomes  exactly  satisfied.    It  can  be  shown  that  under  mild  conditions 

there  exists  a  local  minimum  of  Pq,  x  (p),  for  which 

lim   j:*(p)  =  x  . 

p-=c 

However,  x*{p)  may  not  converge  to  x'  for  any  finite  value  of  p.  As  p  tends  to 
infinity  the  function  becomes  increasingly  more  ill-conditioned.  Moreover,  while, 
for  sufficiently  large  p,  there  exist  local  minima  which  are  near  x  ,  the  penalty 
function  may  be  unbounded  below  for  any  value  of  p.  This  may  cause  the  iterates 
to  move  outside  the  region  of  the  local  minimum. 

A  non-differentiable,  but  well-conditioned,  merit  function  is  the  exact  penalty 
function.  (Han  (1977)) 

Pe(x,p)  ^fix)  +  plk(.v)||i 
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Tht.  advantage  of  Pf  over  Pq  is  that  there  exists  a  finite  value  p*  such  that  for  any 
p>i*,  X  will  be  a  local  minimum  of  Pg.  However,  while  the  existence  of  p  is 
ensi.-red.  the  correct  '/aiue  depends  upon  values  computed  at  x  ,  and  thus  is  not 
knc.  vn  a  priori.  The  algorithm  therefore  may  suffer  from  bad  choices  of  the 
per  :y  parameter.  If  p  is  chosen  too  small,  Pg  may  be  unbounded  below,  or  the 
regii  .\  in  which  x'  is  a  minimum  may  be  very  small.  For  large  values  of  p  the 
pena  /  function  becomes  ill-conditioneci  Since  P^  is  non-differentiable,  standard 
univariate  minimization  methods  designed  for  smooth  problems  cannot  be  used  in 
the  !!•;  search.  Special  algorithms  for  minimszation  of  non-smooth  problems  can  be 
used  lut  these  methods  often  entail  considerable  overhead.  Generally,  a  line 
search  involving  -^  non-differentiable  merit  function  simply  considers  a  sequence  of 
trial  steps,  of  which  the  first  is  unity,  generating  subsequent  trial  points  by  reducing 
the  step  lengths  by  a  constanv  factor  less  than  one,  until  an  acceptable  point  is 
found . 

Since  the  first  order  optimality  conditions  require  x  to  be  a  stationary  point  of 
the  Lagrangian  function  at  .>ii  =  X  ,  it  is  reasonable  to  assume  that  the  Lagrangian 
function  itself  can  be  used  to  measure  progress  towards  a  solution.  However,  while 
,r*  ,s  a  stationary  point  of  the  Lagrangian,  it  is  not  necessarily  a  minimum  The 
seco- d  order  optimality  conditions  specify  that  the  projected  Hessian  of  the 
Lagiiingian  be  at  least  positive  semi-definite.  This  indicates  that  x'  is  a  minimum  of 
the  .agrangian  on  the.  tangent  space  to  the  active  constraints  It  is  possible, 
the. afore,  to  augment  rhe  Lagrangian  function  by  adding  a  term  that  alters  the 
fmK:tion  in  the  range  space  of  A(x')^  so  as  to  cause  x'  to  be  a  local  minimum  Such 
an  lugmented  Lagrangian  penalty  function  (Fletcher  (1981),  Schittkowski  (1981))  is 
giv  n  by 

P.4(x,X,,p)  =  fix)  -  X^c^(x)  +  ^c,.{xVc„{x) 
As  is  the  case  for  the  exact  penalty  function,  there  exists  some  threshold  value,  p. 
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such  that  X  will  be  a  minimum  of  P^  for  all  p>p.  In  addition,  such  merit  functions 
have  the  advantage  of  being  continuously  differentiable.  It  should  be  emphasized 
that  jf  is  a  minimizer  of  the  augmented  Lagrangian  only  when  X  =  X*.  Methods 
incorporating  an  augmented  Lagrangian  merit  function  are  thus  particularly 
sensitive  to  poor  multiplier  estimates. 

The  exact  penalty  merit  function  has  the  advantage  of  not  being  overly 
sensitive  to  multiplier  estimates.  Global  convergence  can  generally  be  proved  for 
methods  which  reduce  an  exact  penalty  function  at  each  iteration.  However, 
insisting  upon  a  decrease  in  an  exact  penalty  function  can  result  in  step  lengths  of 
unity  being  rejected,  even  in  a  very  small  neighborhood  of  the  solution.  This 
phenomenon,  known  as  the  "Maratos  effect",  prevents  some  algorithms  from 
converging  superlinearly  (Mayne  (1979),  Chamberlain  et.  al.  (1982)).  Chamberlain 
et.  al.,  Fletcher  (1982,  1985)  and  Fukushima  (1985)  discuss  techniques  for  avoiding 
the  Maratos  effect. 

Global  convergence  of  methods  using  an  augmented  Lagrangian  penalty 
function  involves  proving  convergence  of  the  multipliers.  Therefore,  Schittkowski 
(1982)  presents  a  globally  convergent  modification  of  the  augmented  Lagrangian 
merit  function  in  which  the  the  line  search  is  performed  over  x  and  X,  rather  than 
over  X  alone.  The  merit  function  used  in  NPSOL  (Gill  et.  al.  (1985c)),  based  on 
that  of  Schittkowski,  is  a  smooth  augmented  Lagrangian  function  of  the  form 

m  1     m 

Pf^ix,K,s,p)  =  fix)  -  2^,(c,  -  i,)  +  72p,(c,-5,)^ 
where  s  are  slack  variables  introduced  in  order  to  keep  the  function  continuously 
differentiable      even      at      constraint      bounds.      The      line      search      evaluates 
P/^(jr*  -i-  ap*",  X*  -I-  a8X*,  s,  p)  to  produce 

The  quantity  8X*  is  defined  by  8X  =  ti*  -  X*^,  where  ti*  are  the  QP  multipliers. 


\'^'     ^     X* 
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Hence,  when  a  step  of  length  one  is  taken,  the  QP  multiplier  estimates  are  accepted. 
In  a  neighborhood  of  the  solution  a  steplength  of  unity  will,  in  fact,  produce  an 
acceptable  point.  Such  a  ment  function  thus  enforces  global  convergence  while 
maintaining  fast  local  convergence  It  should  be  emphasized,  however,  that  this  line 
search  depends  critically  upon  rhe  QP  multipliers  It  is  expected  that  whenever  a 
unit  step  is  laken,  the  QP  multipliers  provide  a  good  estimate  of  the  Lagrange 
multipliers  at  the  new  point. 

Lei  m  denote  some  merit  function.  The  steplength  a*  is  chosen  so  that 
mix''  ^-  a*p*)  is  "sufficiently  lower"  than  m{x''"}.  This  "acceptable"  reduction  in  m  is 
enforced  by  requiring  a''  to  satisfy  the  Goldstein—  Armijo  conditions 

mix^  +  ct^'p*)  <  ot(x*')  +  a'^yVmix'^fp''    -y  €  (0,1)  (3.1a) 

Vm(r*  +  a^'pYp''  >  hVmixYp^  8  €  (-y.l)  (3.1b) 

Intuitively,   the   first  condition   forces   a    sufficient  rate   of  decrease   of  the   merit 

function.  However,  if  arbitrarily  small  steps  are  allowed,  the  method  may  not 
converge.  The  second  condition,  therefore,  excludes  steps  that  are  too  small 
relative  to  the  initial  decrease  of  m  along  the  descent  direction. 

The  fast  local  convergence  rates  of  quasi-Newton  methods  indicate  that  the  full 
step  a*  =  ij  should  be  taken  whenever  possible  The  line  search  is  therefore 
generally  staned  with  a*^  =  1  and  in  ttie  case  that  x*^  +  p''  is  not  acceptable,  some 
"backtracking"  is  done  by  reducing  a'  until  an  acceptable  point  x  +  a  p  is  found 
For  smooth  functions  m,  the  backtracking  can  be  carried  out  by  constructing  a  local 
mode!  to  m(a)  ^  !n(x^  '  ap*)  based  on  the  known  values  of  m(0),  mi  1)  and  m  (0) 
and  choosing  a^  to  minimize  the  mode)  If  the  resulting  step  is  still  unacceptable, 
the  new  information  obtained  about  m  is  used  again  to  construct  another  local  model 
of  m  .\n  easier  means  of  backtracking  and  one  that  is  generally  used  for  non- 
smooth  merit  functior.s  m,  is  simply  to  repeatedly  multiply  a*  by  (,,  for  some  fixed 
C  ^  (0,1)    until   an    acceptable    steplength    is    found      A    backtracking   line    search 
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generally  avoids  excessively  small  steps,  and  therefore  only  condition  (3.1a)  must 
be  considered  in  choosing  a*^.  For  a  more  thorough  discussion  of  step  selection  see 
Dennis  and  Schnabel  (1983). 

Since  our  algorithms  were  developed  by  modifying  the  NPSOL  package,  and 
we  have  tried  to  deviate  from  the  original  code  as  little  as  possible,  we  prefer  to  be 
able  to  use  the  NPSOL  merit  function  Pf^.  As  mentioned  above,  use  of  this  merit 
function  avoids  the  "Maratos  effect"  and  enforces  global  convergence  for  the 
original  NPSOL  algorithm. 

We  would  have  liked  to  show  that  the  global  convergence  associated  with  the 
use  of  the  merit  function  P,v  can  be  extended  to  our  modifications  of  the  original 
method.  However,  as  noted  above,  this  merit  function  depends  critically  upon  the 
QP  multipliers.  Since  our  non-IQP  methods  do  not  solve  a  full  IQP  subproblem,  the 
QP  multipliers  associated  with  a  solution  of  one  of  our  subproblems  are  not  reliable 
estimates  of  the  Lagrange  multipliers  of  the  original  problem.  In  fact,  since  the 
non-IQP  subproblems  do  not  allow  for  deleting  constraints  from  the  working  set  at 
the  start  of  the  QP  (the  PQP  and  EQP  methods  cannot  delete  any  constraints  at  all, 
while  the  MQP  can  delete  only  those  constraints  that  have  been  added  within  the 
current  QP)  the  QP  multipliers  may  have  the  wrong  sign.  In  that  case,  obviously  it 
is  not  logical  to  use  a  merit  function  which  would  set  the  new  multiplier  estimate, 
X**^  to  X.*  -I-  a''(Ti*  -  \.*^).  For  our  IQP  method,  the  QP  multipliers  will  all  be  of 
the  correct  sign.  However,  since  the  multipliers  are  defined  by  the  full  Hessian 
matrix,  (see  (1.9)),  the  QP  multipliers  associated  with  our  IQP  subproblem,  in 
which  the  full  Hessian  is  approximated  by  "unprojecting"  the  projected  Hessian 
approximating  matrix  as  in  (2.12),  may  be  poor  estimates  of  the  Lagrange 
multipliers  of  the  original  problem. 

We  have  considered  various  ways  of  adapting  the  NPSOL  merit  function  for 
use  with  our  modified  methods.    One  possibility  is  to  use  the  least  squares  multiplier 
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estimates  at  x'^  +  p''  in  place  of  r\''  These  multiplier  estimates  may  be  better 
approximations  to  the  Lagrange  multipliers  at  x*^^  (see  Section  3.1i.  However. 
there  is  no  restriction  on  the  signs  of  the  least  squares  multipliers  and  therefore 
they,  too,  may  have  the  wrong  sign  Moreover,  far  from  the  solution,  it  may  be 
necessary  to  take  a  small  step  along  p',  and  the  least  squares  estimates  at  the  full 
step  a:*  +  p*  may  be  totally  unreliable. 

As  a  means  of  avoiding  the  use  of  a  multiplier  estimate  with  the  wrong  sign, 
we  consider  simply  setting  8X*  to  zero  whenever  ti*^  has  the  wrong  sign.  Although 
this  keeps  us  from  revising  the  multiplier  estimates  in  the  wrong  direction,  if  we 
consistently  set  8X*  to  zero  we  will  not  be  modifying  the  multiplier  estimates. 

Near  the  solution  we  expect  that  the  working  set  will  be  an  accurate  predicition 
of  the  active  set  /'  In  that  case,  all  of  our  methods  find  the  solution  to  the 
subproblem  in  one  step  and  the  QP  multiplier  estimates  are  all  of  the  correct  sign. 
In  that  ca-e,  the  NPSOL  merit  function  generally  performs  very  well.  In  fact,  for 
many  problems,  the  working  set  becomes  a  good  prediction  of  the  true  active  set  J 
even  lar  from  the  solution  and  the  NPSOL  merit  function  can  be  used. 

Because  of  the  difficulties  associated  with  using  the  merit  function  Pf^-.  we  also 
include  the  option  of  using  the  exact  penalty  function  P^.  As  noted  above,  this 
merit  function  is  extremely  sensitive  to  the  choice  of  the  penalty  parameter  p,  and 
fast  local  convergence  may  be  inhibited  by  the  "Maratos  effect".  However,  for 
problems  for  which  we  cannot  use  the  NPSOL  merit  function ,  we  have  found  that  in 
practice  the  exact  penalty  function  can  provide  good  results. 

In  later  sections  we  discuss  problems  for  which  a  "natural"  merit  function  is 
available.  Such  problems  are  characterized  the  fact  that  the  objective  function  itself 
can  be  used  to  measure  progress  towards  the  solution.  Since  we  consider  problems 
of  this  nature,  we  include  an  option  in  our  program  to  use  the  objective  function 
itself  as  the  merit  function.  Clearly,  in  these  cases  there  is  no  difficulty  associated 


-54- 

with  the  choice  of  multiplier  estimates  or  penalty  parameter. 

In  summary,  in  implementing  our  algorithms  we  allow  the  user  the  option  of 
choosing  the  merit  function.  The  first  option  is  to  use  the  NPSOL  merit  function 
and  set  8X.|  to  zero  whenever  ti*  has  the  wrong  sign;  the  second  option  is  to  use  an 
exact  penalty  function  and  the  third  option  is  to  use  the  objective  function  itself  as 
the  merit  function.  Numerical  results  which  demonstrate  the  performance  of  these 
options  will  be  presented  in  Chapter  5. 


3.3  Choice  of  a  Basis  for  the  Null  Space  of  A(x). 

Many  techniques  for  constrained  optimization  involve  the  use  of  a  matrix  Zix), 
which  constitutes  a  basis  for  the  null  space  of  the  constraints  at  that  point.  Given 
A(x),  the  fXn  transpose  of  the  Jacobian  of  the  active  constraints,  Z(.r)  is  constructed 
such  that 

A(x)Z(x)  =  0 

Z(xfZix)  =  /„_, 

The  LQ  factorization  of  A(x)   is  a  stable  and  efficient  method  of  forming  Z(x). 

However,  while  A{x)  is  generally  a  continuous  function  of  .v,  Coleman  and  Sorensen 

(1984)   show  that  the  LQ   algorithm  does  not  necessarily  produce  a  continuously 

varying  Zix).    In  fact,  Byrd  and  Schnabel  (1984),  apply  results  from  topology  to 

prove  that  Z  cannot  in  general  be  defined  as  a  continuous  function  of  A(.r),  nor  is  it 

possible  in  general  to  construct  Z  as  a  continuous  function  of  x. 

The  discontinuity  of  the  Z  obtained  by  means  of  the  LQ  factorization  of  A 
arises  from  the  choice  of  sign  associated  with  each  Householder  transformation. 
The  sign  bit  is  chosen  in  such  a  way  as  to  minimize  the  effects  of  cancellation  in  the 
computation  of  the  Householder  matrix.  Coleman  and  Sorensen  present  various 
proposals  for  modifying  the  basic  LQ  method  so  as  not  to  involve  changes  in  the 
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sign  )f  the  transformations.  The  first  method  involves  reordering  the  columns  of  A 
so  as  to  keep  the  sign  bit  the  same  for  all  x  within  some  neighborhood  of  .v  .  The 
secord  proposal  is  motivated  by  the  observation  (due  to  Parlett  (1971))  that  if  the 
comp  ration  of  the  transformation  is  computed  in  an  appropriate  manner,  the 
necessity  of  using  the  sign  bit  as  a  safeguard  is  obviated  Thirdly,  Coleman  and 
Scren.-en  suggest  the  use  of  elementary  rotation  matrices  rather  than  the  standard 
elementary  reflectors  in  forming  the  factorization 

Wnile  these  methods  are  convergent  in  some  neighborhood  about  x  ,  the 
resuhi  of  Byrd  and  Schnabel  indicate  that  for  a  large  class  of  matrices,  Z{^)  must 
have  -a  discontinuity  in  any  ball  around  A  of  radius  larger  than  its  smallest  singular 
value  Thus  the  resuUs  of  Coleman  and  Sorensen  may  be  applicable  only  in  a  very 
small  neighborhood  of  v 

The  significance  of  these  results  lies  in  the  implications  they  have  for  methods 
which  maintain  an  approximation  to  the  projected  Hessian  of  the  Lagrangian 
function.  The  projected  Hessian,  Z*'W(;c*,\*)Z'',  is  dependent  upon  the  choice  of 
basis  for  the  null  space.  The  fact  that  in  general  Z*  cannot  be  constructed  as  a 
continuously  varying  function  of  v*  would  make  it  difficult  to  prove  convergence 
results  for  such  methods. 

Byrd  and  Schnabel  suggest  modifying  the  standard  quasi-Newton  updating 
algor' iims  in  such  a  way  as  to  ensure  that  the  iterates  are  independent  of  the  choice 
of  Z'  Their  proposal  is  to  regard  the  recurred  matrix,  S*,  generally  viewed  as  an 
appn  .imation  to  the  projected  Hessian,  as  implicitly  providing  information  about 
the  en:ire  Hessian.  Given  S*,  an  approximation  to  Z*'vV(.r*,\*)Z*,  "unproject"  S*  to 
form 

M*  =  Z'^B^Z'"'  +  p*(/-Z*Z*') 
(wher.    p*  is  an  estimate  of  ||W(x*,X.*)  |!  ).    Now  that  the  full  Hessian  has  been 

retrieved,  the  quasi-Newton  update  can  be  applied  to  the  full  Hessian  projected  onto 
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the  subspace  defined  by  Z''*^.    Form 

z'^'''m'z'*'  =  Z'^^'Z^B^Z^'Z'*'  +  p*(/-z*"i'z*z*'z*^i) 

and  apply  the  quasi-Newton  update  in  the  projected  space  to  this  matrix  to  form 
S*"^^  (Note  that  when  Z'^=Z'''^^,  this  formulation  reduces  to  the  original  methods.) 
Although  this  suggestion  does  produce  iterates  that  are  independent  of  the  definition 
of  the  null  space  basis,  it  should  be  noted  that  no  convergence  results  have  been 
proved  as  yet  for  algorithms  which  incorporate  a  subspace  shift. 

A  different  approach  is  adopted  by  Gill,  Murray,  Saunders  and  Wright(  1983d). 
They  present  a  method  in  which  Z  is  computed  by  updating  the  TQ  factorization  of 
the  previous  iteration.  This  Z  is  continuous  in  the  neighborhood  of  a  point  where  A 
has  full  rank;  if  Z  is  computed  at  a  sequence  of  points  converging  to  x  ,  the 
sequence  {Z*}  converges  as  well.  However,  the  same  results  do  not  apply  to  the  full 
matrix  Q.  Moreover,  it  should  be  noted  that  this  Z*  is  generated  as  a  continuous 
function  of  A*  and  the  TQ  factorization  of  a''~^,  but  that  Z*  is  not  a  function  of  A* 
nor  of  ,r*. 

In  a  later  paper  (1985),  Gill  et  al.  extend  these  results  and  derive  explicit 
bounds  on  the  change  in  Z  with  respect  to  perturbations  in  A.  They  introduce  a 
class  of  regularized  Householder  transformations,  which  differ  from  the  standard 
Householder  transformations  by  a  reversal  of  the  sign  of  one  row.  By  using 
regularized  Householder  transformations  to  update  the  TQ  factorization  in  place  of 
the  standard  updating  method,  they  are  able  to  derive  a  bound  on  the  change  in  Q. 
Continuity  of  the  full  matrix  Q  can  be  proven  in  exactly  the  same  manner  as  for  Z. 

Since  there  are  no  convergence  results  yet  for  the  updating  method  due  to  Byrd 
and  Schnabel,  and  in  order  to  preserve  the  similarities  between  our  algorithms  and 
NPSOL,  we  prefer  to  use  the  TQ  factorization  update  developed  at  Stanford.  We 
have,  however,  experimented  with  the  subspace  shift  suggestion  of  Byrd  and 
Schnabel  on  a  small  test  problem  set.    Our  results  did  not  indicate  any  significant 
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difference  in  the  performance  of  these  two  methods.  However,  a  more  thorough 
c^  mparison  of  the  two  methods  is  necessary  before  any  conclusion  can  be  reached 
a:'    o  their  relative  worth. 


3,4  Ouasi-Newton  Updates 

Given  S*,  a  positive  definite  approximation  to  Z*^'W*Z*,  we  wish  to  apply  a 
stand,  d  updating  method  sucli  as  DFP  or  BFGS  to  obtain  S**^  For  the 
ancon-:ramed  case,  the  BFGS  update  is  defined  by 


and  the  DFP  update  is  given  by 


5^-1  =  g/:  +  JUL-  _   "  ^  ^    "  (3.2) 


where 


Empidcal  results  and  the  analysis  of  Powell  (1986)  and  Byrd  and  Noceda!  (1986) 
mdicate  that  the  BFGS  update  is  preferable 

When  the  updates  are  to  be  performed  in  the  projected  space,  alternative 
de  nitions  must  be  used  for  .5*  and  >*  Nocedal  and  Overton  M985),  in  discussing 
t'n    equality  constrained  problem,  suggest  the  following  choices; 

s'  =  Z''^'\x'''^  -  x'')  (a) 

or    s'  =  Z*Yr*^l  -  x')  (b) 

/  =  Z*-i'-(/-i-  ,/-  A^^X^i))  (c) 

or   y'  =  Z'-'\g'''  -  ^g'-  A'\i))  .                                       (d) 

or   y'  =  Z'\g'^'-  A'-'\i^'-  g')  (e) 
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or   y*  =  Z*-'^/-l  -  Z*V*  (f) 

The  definitions  for  y*  are  derived  from  considering  the  change  in  the  gradient  of  the 

Lagrangian  function  (1.7).    For  example,  (c)  is  equivalent  to 

For  the  inequality  constrained  problem,  the  number  of  active  constraints  may 
change  from  one  iteration  to  the  next.  Thus  Z*  and  Z*"^^  and  hence  fi*  and  S*"^^ 
may  not  be  of  the  same  dimension.  In  order  to  extend  these  definitions  of  .?*  and  y* 
to  the  more  general  problem,  we  replace  Z*  in  (a)-(f)  by  Z^  which  is  defined  as  a 
basis  for  the  null  space  of  A^,  i.e.,  the  matrix  whose  rows  consist  of  the  gradients  of 
the  constraints  that  are  active  at  a:*"^^  but  evaluated  at  .r*.  If  more  constraints  are 
active  at  x'''^^  than  at  .r*,  the  transformations  applied  to  update  the  TQ  factorization 
are  applied  to  the  Cholesky  factor  of  S*  as  well  (see  Section  2.2).  When  fewer 
constraints  are  active  at  jf*"^\  the  approximation  to  the  projected  Hessian  matrix 
must  be  expanded  by  adding  a  row  and  a  column.  When  second  derivative 
information  is  available  for  the  full  space,  the  new  column  can  be  computed  by 
projecting  this  information  into  the  new  space.  For  our  purposes,  we  simply  use  a 
positive  multiple  of  the  identity  matrix  in  place  of  the  full  Hessian. 

Unlike  unconstrained  optimization,  there  is  no  guarantee  that  5*  y*  >  0. 
Therefore,  Womersley  and  Fletcher  ( 1984a, b)  use  a  modification  of  the  BFGS 
method  proposed  by  Powell  (1978a)  which  guarantees  hereditary  positive 
definiteness.   Powell's  modification  consists  of  replacing  y*  in  (3.2)  by  the  vector 

r*  =  ey*  +  (1  -  8)fi*5*, 
where  0  is  defined  by 


e 


1,  /yso.25*'sv. 


5*  y*  <  0.25*  fi*j*. 


This  modification  to  the  update  prevents  the  determinant  of  B        from  being  less 
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:han  0.2  of  the  determinant  of  B*.  (The  constant  0.2  was  chosen  as  a  reasonable 
size,  but  not  because  of  any  theoretical  properties.) 

However,  using  the  above  definitions  of  5*  and  >•*,  it  is  possible  to  construct 
examples  where  ||  .?*  ||  is  small  while  \\  y''  \\  is  not.  The  Powell  update  does  not 
avo'd  this  problem,  since  it  may  be  the  case  that  ^^'y*  >>  .J*'i5*5*.  In  this  event,  the 
updates  may  be  undefined  or  may  cause  the  Hessian  approximation  to  blow  up 
Coleman  and  Conn  (198^)  circumvent  this  problem  by  introducing  an  intermediate 
point,  X*  =  j:*  +  Z'^p' .  and  defining 

This   pncedure  ensures  that   \\  s''  \'.  cannot  be  too  small  compared  to   \\y^  ||,  but 
requires  an  additional  set  of  gradient  evaluations  at  each  iteration. 

.V.  cc,  al  and  Overton  suggest  an  alternative  strategy.      If  i|5*  ||  is  too  small, 
they  simph  do  not  perform  the  update     Let 

v*=  Y'*'\x'"  -  x') 
or   V*  =  K' U*"^  -  x") 
corresponding  to  whether  (a)  or  (b)  is  used  to  define  5*.    The  update  is  done  if  and 

only  'f 


■*here  w  and  v  are  appropriately  chosen  constants.    Under  certain  conditions  this 
update  rule  will  not  impede  local  two-step  superlinear  convergence 

In  practice,  we  have  found  that  it  is  generally  preferable  not  to  skip  updates 
The  update  rule  of  Nocedal  and  Overton  may  be  too  restrictive  and  may  cause  the 
algorithm  to  skip  updates  even  in  a  situation  in  which  updating  the  projected 
Hessian  would  not  cause  dny  difficulty. 

Ideally,  we  would  like  to  use  some  update  rule  which  would  guarantee  that  the 
eigenvalues  of  the  projected  Hessian  approximating  matrix  lie  within  given  bounds. 
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Calamai  and  More  (1985)  have  recently  developed  a  quasi-Newton  method  for 
solving  systems  of  nonlinear  equations,  based  on  the  Broyden  update,  which 
preserves  bounds  on  the  elements  of  the  Jacobian  matrix.  However,  it  is  not  clear 
how  to  extend  this  method  to  preserve  positive  definiteness  nor  how  to  apply  it  in 
the  projected  space. 

Instead,  we  consider  another  simple  updating  rule  based  on  an  idea  of  Sargent 
(private  communication).  We  use  the  Powell  modification  of  the  BFGS  update  if 
this  preserves  some  bounds  on  estimates  of  the  eigenvalues  of  the  projected  Hessian 
matrix.  Since  we  recur  the  Cholesky  factor  of  the  projected  Hessian  matrix,  /?*,  we 


use  the  values  r^^^,  r^,„  and 


^2 
'max 


as  estimates  of  the  maximum  and  minimum 


eigenvalues  and  the  condition  number  of  the  projected  Hessian,  respectively,  where 
''max  =  rnax  ||/?*(/,/)||  and  r^^^  =  min  ||/?*(/,/) ||.    If  the  estimates  of  the  eigenvalues 

of  the  new  projected  Hessian  approximation  do  not  lie  within  the  set  bounds,  we  do 
not  update.  This  update  rule  is  less  restrictive  than  the  Nocedal  and  Overton  rule 
and  yet  will  not  allow  the  approximating  matrix  to  blow  up. 

Nocedal  and  Overton  discuss  an  example  for  which  ||5*||  «  ||y*||  and 
performing  the  update  would  cause  the  norm  of  the  reduced  Hessian  approximation 
to  become  very  large,  while  use  of  their  updating  rule  would  avoid  this  difficulty. 
We  note  that  our  updating  rule  would  provide  the  same  results.  Another  example 
presented  by  Nocedal  and  Overton  demonstrates  the  sensitivity  of  their  update  rule 
to  the  proper  choice  of  co.  An  incorrect  choice  of  o)  can  allow  the  update  to  be 
performed  and  cause  the  projected  Hessian  approximation  to  become  nearly 
singular.  The  advantage  of  our  updating  rule  is  that  any  such  update  which  would 
result  in  a  near-singular  matrix  would  be  skipped. 
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3,5  Dropping  Criteria 

Many  criteria  for  deciding  when  to  drop  a  constraint  from  the  working  set  have 
been  studied.  On  one  extreme,  exact  minimization  on  the  current  subspace  can  be 
required  before  a  constraint  is  dropped.  On  the  other  extreme,  some  criteria  call 
fo  deleting  a  constraint  whenever  its  associated  Lagrange  multiplier  estimate  is 
negative.  In  general,  requiring  subspace  minimization  tends  to  cause  the  algorithm 
to  take  a  greater  number  of  iterations,  since  minimizing  on  the  wrong  subspace  does 
■not  necessarily  produce  progress  toward  the  solution.  However,  an  algorithm  which 
may  drop  constraints  at  every  step  may  suffer  from  zigzagging  An  example  due  to 
Wolfe  ri972)  shows  that  such  a  method  may  converge  to  a  non-optimal  point. 

When  there  is  more  than  one  constraint  with  a  multiplier  estimate  of  the  wrong 
sign,  ;t  is  usual  to  delete  the  constraint  corresponding  to  the  largest  multiplier 
estimate  in  absolute  value  with  an  incorrect  sign  (see  Fletcher  (1981)).  However, 
we  have  found  that  where  the  constraints  are  nonlinear,  this  may  not  be  the  best 
choice  A  nonlinear  constraint  may  be  included  in  the  working  set  and  yet  may  not 
oe  exactly  satisfied.  In  fact,  the  linearized  constraint  may  be  exactly  satisfied  while 
the  current  iterate  is  actually  infeasible  with  respect  to  the  original  nonlinear 
co^i-^raint. 

Another  related  issue  is  how  much  to  modify  the  working  set  at  each  iteration. 
Allowing  only  one  constraint  to  be  deleted  at  each  iteration  may  cause  the  method 
to  take  man>  ^teps  until  the  working  set  is  correctly  determined;  dropping  many 
constraints  too  early  may  cause  the  method  to  take  extra  steps  in  order  to  pick  them 
up  again.  In  practice,  we  have  found  that  constraints  can  be  deleted  from  the 
working  set  until  the  gradient  projected  onto  the  new  space  is  no  longer  small. 

In  later  sections  we  will  discuss  the  theoretical  properties  associated  with  an 
additional  condition  on  dropping  rules  proposed  by  Byrd  and  Shultz  (1984).  This 
condition  is  that  no  constraint  may  be  deleted  from  the  working  set  if  the  previous 
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iteration  involved  adding  a  constraint,  unless  the  current  point  is  a  non-optimal 
vertex.  Enforcing  this  additional  condition  allows  us  to  apply  some  of  the 
theoretical  results  of  Byrd  and  Shultz  to  our  methods.  However,  in  practice  the 
additional  restriction  on  dropping  constraints  may  cause  the  algorithm  to  take  extra 
steps. 

In  our  implementation,  we  allow  the  user  the  option  of  enforcing  the  Byrd  and 
Shultz  rule,  or  of  disregarding  it,  by  setting  an  appropriate  flag.  The  user  may  also 
control  how  far  to  minimize  on  the  current  subspace  before  allowing  a  delete 
operation  by  setting  the  tolerance  level.  If  desired,  the  user  may  allow  a  default 
tolerance  level  to  be  used,  which  is  computed  by  the  algorithm  and  is  dependent  on 
the  scaling  of  the  problem.  If  the  projected  gradient  is  less  than  the  tolerance  level, 
the  method  may  delete  one  or  more  constraints  until  the  gradient  projected  onto  the 
larger  space  is  above  the  tolerance  level.  The  constraints  to  be  deleted  are  chosen 
to  correspond  to  the  largest  least  squares  multiplier  estimates  in  absolute  value  with 
an  incorrect  sign,  chosen  from  among  those  constraints  that  are  satisfied. 

Because  of  the  difficulties  involved  in  choosing  an  appropriate  merit  function 
(see  Section  3.2),  it  may  be  the  case  that  in  order  to  find  an  acceptable  point  for  the 
merit  function  along  the  search  direction  /)*  the  step-length  must  be  extremely  small. 
If  this  occurs,  unless  the  working  set  is  modified,  the  new  iterate,  jr*^^  will  be  very 
close  to  x''  and  p'^*^  will  be  close  to  p''  and  the  method  will  keep  taking  extremely 
small  steps  until  it  moves  far  enough  away  from  x'' .  Therefore,  if  we  find  that  a 
very  small  step  is  taken,  we  immediately  try  to  modify  the  working  set  by  deleting  a 
constraint,  regardless  of  the  size  of  the  projected  gradient. 
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Chapter  4 
Convergence  Results  for  a  Restricted  Class  of  Problems 


4.1  General  Nonlinear  Programming  Problems 

For  the  purposes  of  this  chapter,  we  consider  the  problem  to  be  again  in  the 
form  (NP),  i.e.  with  all  constraints  written  with  lower  bounds. 

Global  convergence  of  SQP  methods  for  the  general  inequality  constrained 
nonlinear  programming  problem  has  been  shown  by  Han  (1977),  Schittkowski 
(1982)  and  Gill  et.  al.  (1986,  to  appear)  for  IQP-based  methods,  and  by  Coleman 
and  Conn  (1982b)  for  an  EQP-based  method.  Schittkowski  and  Gill  et.  al.  base 
their  proofs  of  global  convergence  upon  the  use  of  an  augmented  Lagrangian  merit 
function  where  the  line  search  is  performed  over  both  x  and  X.  (see  Section  3.2). 
However,  as  discussed  in  Section  (3.2),  such  a  merit  function  cannot  always  be  used 
for  our  methods.  Global  convergence  of  Coleman  and  Conn's  method  is  predicated 
upon  an  extra  evaluation  of  the  objective  and  constraint  functions  at  each  step. 

Han's  proof  most  readily  can  be  extended  to  our  IQP  method.  If  we  assume, 
as  does  Han,  that  S*  has  bounded  eigenvalues  for  all  k,  and  that  the  parameter,  p, 
of  the  exact  penalty  merit  function  is  chosen  correctly,  then  we  can  apply  Han's 
Theorem  (3.1)  to  show  that  p'' ,  the  direction  produced  by  our  IQP  method,  is  a 
descent  direction  for  the  exact  penalty  merit  function.  To  ensure  "sufficient 
decrease"  in  the  merit  function  Han  chooses  the  steplength  a*  so  as  to  satisfy 

f£(x*  +  aV*,P)  ^    min    Pgix''  +  ap*,p)  +  €*, 

Osa£l 

X 

where  {e*}  is  a  sequence  of  nonnegative  numbers  satisfying    2)  €*  <  o=.    If  we  use 

i  =  0 
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Han"?  condition  for  the  line  search  and  make  the  additional  assumption  that  strict 
complementarity  holds  at  x' ,  we  can  apply  Han's  Theorem  (3.3)  to  show  that  our 
IQP  method  is  globally  convergent  for  problems  in  which  the  constraint  functions 
are  con  ave. 

Usually,  proving  global  convergence  for  the  general  problem  is  difficult  and 
deperid>  upon  the  choice  of  an  appropriate  merit  function.  As  discussed  in  Section 
3.2,  tht  question  of  how  to  select  an  appropriate  merit  function  is  itself  a  difficult 
problem  Therefore,  in  ihe  following  sections  we  will  restrict  our  attention  to 
certain  classes  of  problems  for  which  this  difficulty  does  not  arise. 


4.2  Minimum  1^  and  li  Problems 

For  general  nonlinear  programming  problems,  choice  of  an  appropriate  merit 
function  is  difficulr  However,  there  is  a  class  of  problems  for  which  a  "natural" 
merit  function  is  available  Let  f{x):  R''-R'"  be  a  twice  continuously  differentiable 
function.   Consider  the  problems 

min    i(x]  =  \\f(x)  111  (l^) 

xiR" 

and 

mm   y(x)^     \\fix)  |U  ('/) 

xiR" 

where 

II /w  iii=  ii/,w  I 

and 

Wfix)  \L  =    max  !/,(x)  !. 
These  problems  arise  from  data-fitting  applications,  where  one  attempts  to  fit  the  m 
observations  with  a  nonlinear  function  of  n  variables.    The  problem  consists  of 
selecting  x  so  that  the  fit  is  optimal  in  some  sense.  By  minimizing  the  /^  norm  x  is 
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chosen  so  that  the  average  deviation  of  the  model  from  the  data  is  minimized,  while 
the  infinity  norm  is  used  to  minimize  the  maximum  deviation. 

Problem  (l^P),  also  known  as  the  minimax  problem,  is  equivalent  to  the 
following  nonlinear  programming  problem  in  n  +  1  variables: 

min         8  (PJ 

(;t6R",  8€Ri) 

s.t.    S  -  f,{x)  >  0 

8  4-  fiix)  >  0,  J  =  l,...,m 

Similarly,  by  introducing  m  extra  variables,  (liP)  can  be  transformed  into 

min  2 "/  (^/) 

S.t.    u,  -  Mx)  >  0 

u,  +  f,(x)  >  0,  /•  =  l,...,m 

Any  algorithm  for  solving  nonlinearly  constrained  optimization  problems  can 

be  used  to  compute  the  solution  to  these  formulations  of  (/jP)  and  {l:^P).    These 

problems  have  the  additional  property  that  feasibility  can  always  be  maintained. 

Specifically,  for  any  value  of  .v,  8  and  u  can  be  chosen  by 

8  =  II /(;c)  II. 
"/  ~  I  fi^^)  I    /  =  1,.. .,/" 
The    constraint    functions,    therefore,    can    be    ignored    when    measuring    progress 

towards  the  solution,  and  the  objective  functions  themselves  can  be  used  directly  as 
the  merit  function. 

We  feel  that  this  class  of  problems  provides  a  setting  in  which  the  performance 
of  our  algorithms  can  be  measured.  The  existence  of  a  natural  merit  function  allows 
us  to  focus  our  attention  on  the  underlying  methods  without  concerning  ourselves 
with  anomalies  introduced  by  a  poor  choice  of  merit  function.  The  viability  of  any 
method  for  solving  inequality-constrained  nonlinear  programming  problems  depends 
upon  its  ability  to  correctly  identify  the  active  constraints.  Finding  the  solution  to 
minimax  and  minimum  /]  problems  typically  involves  many  changes  in  the  working 
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set.    Performance  of  our  algorithms  on  problems  of  this  nature  will  provide  an 
indication  of  their  capabilities  with  regard  to  changes  in  the  working  set. 

Following  the  standard  definition  of  the  QP  subproblem,  the  IQP  formulation 
of  the  subproblem  for  (Pj,}  would  be  given  by 

min     y,p'^Hp  +  el^^  (4.1) 

s.t.      A^p  ^  -c 
where    //   is   an   (n  +  l)X(n-l- 1)    positive   definite   matrix    which   approximates   the 

Hessian  of  the  nonlinear  problem  and  A  and  c  are  defined  appropriately.    The  EQP 

formulation  follows  simikrly.    A  standard  SQP  method  would  then  search  along  the 

direction  p  and  set  the  new  iterate  to  t*"*"^  =  x'^  +  a.p .  However,  such  a  method 

ignores    the    knowledge    of  the    special    structure    of   the    objective    function    and 

constraints 

A  numbei  of  algorithms  for  solving  the  nonlinear  minimax  problem  which 
exploit  the  special  structure  of  the  problem  have  been  studied  (Conn(1979),  Murray 
and  Overtone  1980)  and  Overton  (1982)),  While  these  methods  are  more  efficient 
for  this  problem,  they  do  not  necessarily  conform  to  the  general  framework  for 
nonlinear  programming  problems.  We  are  interested  in  studying  minimax  problems 
as  a  special  case  of  the  general  nonlinear  programming  problem,  and  are  therefore 
mterested  only  m  such  methods  as  retain  the  basic  structure  of  the  general  problem, 
ye:  allow  for  some  slight  modification  in  view  of  the  special  structure  of  the 
minimax  problem 

We  consider  two  slight  variations  of  the  basic  method  which  exploit  this 
knowledge,  yet  fit  into  the  standard  framework  for  general  nonlinear  programming 
problems.  In  the  first  case,  we  pose  the  QP  subproblem  as  (4.1),  but  ignore  the 
value  of  p„^.l.  Thus,  the  QP  is  posed  in  R"*",  yet  the  line  search  is  done  in  R". 
This  approach  is  motivated  by  the  fact  that  the  minimax  function  depends  only  on 
the  first  n  variables.    The  a  +  lst  variable  is  just  an  approximation  to  the  objective 
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function.  Therefore  the  line  search  should  be  done  in  R"  to  determine  new  values 
for  the  first  n  variables,  and  S'^^^  should  be  defined  by  the  value  of  the  minimax 
function  at  the  new  point. 

For  the  second  variant,  we  note  that  since  the  second  derivative  of  8  is  zero, 
the  Hessian  approximation  in  the  QP  should  be  of  the  form 


H  = 


H    0 
0     0 

Although  H  is  now  indefinite,  H  is  positive  definite  and  Z  HZ  will  still  be  positive 


definite  at  the  solution.    This  second  approach  is  the  one  used  by  Han  (1978,  1981). 
He  defines  the  IQP  subproblem  as 

min        \p'^Hp  +  8  (4.2) 

s.t.     |/,(.0  +  f[{x)p  I  <  8 
From  the  definitions  of  A  and  c.  the  constraint  condition  is  the  same  as  A^p  ^  -c. 

Again,  we  note  that  8  is  not  used  in  the  line  search.    In  the  following  discussion,  we 

will  assume  that  the  subproblem  is  defined  as  (4.2)  or  the  equivalent  formulation 

for  the  EQP,  PQP  or  MQP.    We  stress  that  for  our  methods  the  matrix  H  is  derived 

from  "unprojecting"  an  approximation  to  the  projected  Hessian  (see  Section  1.5). 

Since  our  IQP  method  corresponds  to  that  of  Han,  we  are  able  to  apply  Han's 
results  to  prove  the  global  convergence  of  our  IQP  method  for  the  minimax 
problem. 

The  following  theorem,  due  to  Han  (1981,  Theorem  3.1),  can  be  directly 
applied  to  our  method: 

Theorem  4.1:  If  p*  =  0  then  7(.v'^)  =  S*"  and  the  point  .v*  is  a  stationary  point 
of  Problem  (P^)  with  the  vector  ti*  (defined  as  in  2.9)  as  its  associated  Lagrange 
multiplier. 

Han  then  proceeds  to  show  that  the  search  direction  p*  is  indeed  a  descent 
direction  for  -y  when  it  is  nonzero  (Lemmas  3.2  and  3.3  and  Theorem  3.4  of  Han 
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(1981)).  Again,  these  proofs  hold  for  our  method  as  well.  Han's  Theorem  3.6 
establishes  that  if  p  #0  and  the  matrix  W*  is  positive  definite,  the  steplength  using 
the  minimax  function  y  as  the  merit  function  will  be  some  positive  value  and  can  be 
found  by  a  backtracking  line  search  in  a  finite  number  of  steps. 

These  results  are  used  ".o  prove  the  main  convergence  theorem  (Theorem  3.8  of 
Han): 

Theorem  4.2:  If  yix)  is  bounded  below,  the  Hessians  of /,(.v)  are  all  bounded, 
an.i  if  there  exist  Kj,  sc2  >  0  such  that  Kix'^x  >  x^fi''  x  ^  i<2^^'''  fo""  ^'^  v^R"  and 
every  k,  then  either  the  sequence  {x'^}  generated  by  the  IQP  method  terminates  at  a 
stationary  point  of  Problem  (P_J  or 

lim  m/||p*||  =  0. 


4.3  Linearly  Constrained  Problems 

For  nonlinear  problems  in  which  all  the  constraints  are  linear,  the  problems 
associated  with  choosing  an  appropriate  merit  function  disappear  Since  feasibility 
can  easily  be  maintained  for  linear  constraints,  the  line  search  can  depend  on  the 
objective  function  alone  In  the  following  we  assume  that  all  the  constraints  are 
linear  and  that  the  initial  iterate,  r°,  is  feasible.  Under  these  conditions,  we  show 
that  ail  four  of  our  methods  produce  descent  directions  for  the  objective  function. 
We  shall  also  discuss  an  active  set  strategy  which  guarantees  non-zigzagging  and 
global  convergence  for  our  EQP.  PQP  and  MQP  methods.  Note  that  since  an  IQP 
method  does  not  entail  modifying  the  working  set  during  the  major  iterations,  the 
framework  of  this  convergence  oroof  does  not  apply  to  such  methods. 

Theorem  4. .3:  Let  c'^  =  0.  Then  p",  the  solution  to  the  EQP,  is  a  descent 
direction  for/,  provided  it  is  not  zero. 
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Proof:  Note  that  since  c*  =  0,  no  step  is  taken  to  the  constraints,  i.e.,  Py  =  0. 
The  proof  then  follows  immediately  from  the  definition  of  p*  and  positive 
definiteness  of  fi*: 

^*p*  =  ,?*'^zVr  since  p*  =  0 

=  -g*'z*S*"'z*'^e;>    from    (2.14) 

=  -g'''z''B''"z'''g''      from  (2.8)  and  the  fact  that  p*  =  0 

<  0. 
Recall  that  g''  and  6*  denote  the  gradient  of  the  objective  function  and  the 
projected  Hessian  matrix  at  the  start  of  the  ^th  iteration,  respectively.  The 
approximation  to  the  full  Hessian  matrix,  //*,  is  obtained  by  "unprojecting"  6*  as  in 
(2.6).  As  defined  in  Chapter  2,  the  notation  A*  refers  to  the  matrix  of  transposed 
active  constraint  gradients,  Zf  gives  a  basis  for  the  null  space  of  A,  and  iB,  denotes 
the  projected  Hessian  approximation  at  the  ith  minor  iteration  of  the  ^th  QP 
subprobiem.  Other  quantities  are  defined  similarly.  Note  that  although  the 
approximation  to  the  projected  Hessian  matrix  is  modified  within  the  PQP,  //*  is 
always  defined  in  terms  of  Sq  ^^^  '^  is  always  true  that 

fif  =  Zt'H^zt  (4.3) 

(see  Step  7  of  the  algorithm).    We  assume  that  unless  Sf  is  the  empty  matrix,  it 

satisfies  Kj.v^a:  ^  x'^B\x  ^  k^.v^.v,  for  some  nonzero  Kj  ,K2  and  all  z,  k.  It  follows 
from  (4.3)  that  the  same  is  true  for  H^  when  D  is  chosen  to  be  a  positive  definite 
matrix.  When  it  is  clear  that  all  quantities  are  being  evaluated  at  major  iteration  k, 
we  will  omit  the  superscript. 

Consider  Z,  +  i  obtained  when  a  constraint  is  added  to  the  working  set.  A  row 
is  added  to  A,  to  form  A,^i.  The  TQ  factorization  of  A,  can  be  updated  to  produce  a 
factorization  of  A.  +  j  by  applying  a  sequence  of  plane  rotations  to  Q^.  Since,  in  the 
context  about  to  be  used,  only  the  range  of  Z,  is  important,  we  can  assume,  without 
loss   of  generality,   that   Z,  +  i  =  [  Z  ]    and   Zj  =  [Z  z  ],    where   Z^z  =  0,    although 
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computationally  thi^  may  not  be  the  case. 

Lemma  4.4:  Given  A,,  if  A,*;  is  obtained  by  adding  a  constraint  to  the  working 
et,     and     Z,^i     is     formed     by     updating     the     TQ     factorization     of    A,,     then 

Proof:  The  matrix  product  Z,Zf  is  a  projector  matrix  operating  in  the  rows  of 
ZT^  .  or  equivalently,  on  the  columns  of  Z,_i.  Since  the  columns  of /;^i  already  are 
in  ti.e  range  of  Z,.  the  projection  does  not  change  them. 

A  similar  argument  shows  that  Zj+iZ^Zf-  Zf+i  where  Z,  +  /  is  obtamed  from  the 
upti-u.d  factorization  of  A,  when  /  constraint.s  have  been  added  to  the  active  set  with 
no  rnervening  deletions. 

Lemma  4.5:  Let  Zj  and  Z.^i  be  defined  as  above.  The  matrices  Y,  and  y,^i  are 
similarly  obtained  from  the  TQ  factorizations  of  A,-  and  A,  +  i  as  in  (2.4).  The  matrix 
product  Zj^iY,  =  0. 

Proof:  The  proof  follows  immediately  since  the  columns  of  Z,^;  lie  in  the  range 
of  Z,  and  zJy,  =  0  by  definition.  Similarly,  Z,>;K,-  =  0,  where  Z^  .  is  defined  as 
above 

Theorem  4.6:  .Assume  that  x'^  is  feasible,  c*^  =  0.  and  A*  has  full  rank.  Then 
/?*,  the  solution  of  the  PQP.  is  a  descent  direction  for/  provided  it  is  .lot  zero 

Proof:  Let  us  first  consider  H^x,  the  effect  of  multiplying  a  vector  by  the 
estimated  full  Hessian.    Since 


we  have 


B'     0 

U      D 


H''x  =  Z^B''Z^\  +  Y^DY^\.  (4.4) 

We  will  consider  the  cases  uhere  the  diagonal  block  is  either  zero  or  the  identity 

matrix.    In  the  former  case,  the  second  term  drops  out.    Note  that  in  any  case, 
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within  the  PQP  we  use  the  vector  //*.r  only  when  premultiplying  it  by  Z^  . 
Therefore,  the  second  term  involves  a  matrix  product  of  the  form  Zf  K*  which  must 
be  zero  due  to  Lemma  4.5. 

Since  c*  =  0,  A*  has  full  ranlc  and  the  columns  of  K*  span  the  range  space  of 
/\*  ,  we  have  p*  =  0.    Hence  p  =  2^f/'.'*-    Recall  that  a*  denotes  the  steplength 

taken  at  Hh  minor  iteration  of  the  ^th  QP  subproblem.  Note  that  0  ^  af  <1,  and 
that  if  ctf  =  1,  Zf+igQp  =  0,  so  that  the  subproblem  is  terminated.  In  the  following 
we  will  drop  the  superscript  k.  Consider  p_,: 

gQP,  =  g  +  fiYiPy   from  (2.8) 

=  g  since  p*  =  0  (4.5) 

Bip^   =  ~ZigQp         from  Step  7.11  of  the  PQP  algorithm 

Sqp   ~  Sqp  "•"  ^ifiZip,    from  (2.8)  and  the  definition  of  p 
^iSqp,  =  ZligQP^  +  a-^HZ^,) 

=  Zlg  -  a,ZlHZ,B^'z\g 
=  Zlg  -  a^zlHZi(Z\HZ{)~^z\g  from  (4.3) 
=  2lg- 
^{Zl{Z{Z\  +  Y^Y\)HaxZ\  +  Y{i\)Z^{Z\HZ{)-^Z\g 
since  [Zj  KJ  is  orthogonal 
=  Z\g  -  ■Si{ZlZ{Z\HZ{Z\Z^{Z\HZ^)-^Z\g 
since  Y\Zi  and  Y\Z2  =  0  by  the 
definition  of  Zj  and  Lemma  4.5. 
=  Z\g  -  a^ZlZ.Zlg 
=  Zlg  -  a.{Zlg 
=  {\-a^)Zlg  (4.6) 

=  -(l-ai)Z[^ 
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ZlgQp,  =  Zhgp,  -  ^iZlHZ.iZlHZ.J-'zlgQp^ 
=  (1  -  a.2)ZlgQp_  as  abo%'e. 


Now,  from  the  relationship  between  Zi  si^d  Z3,  we  have 
ZlgQP,  = 


zhop, 

-3.?QP, 


and 


ZlgQP. 


Z\gQP, 

T 
^38QP, 


From  i'^.6) 


Therefore, 


Similarly, 


Z28QP,  =  ( 1   -  (^\)Z2gQP^. 
ZlgQP,  =   (1    -   ^l)ZlgQP^ 

ZlgQP^=  (1  -  ai)(l  -  K2)Zh- 

Sa/'.',  =  -zlgQP, 

=  -tl-a,)(l-a2)Zlg. 

Bj>z  =  -(  nn-a,)  )Zfg. 


We  arc  now  able  to  show  that  /?*  is  a  descent  direction: 

~    ~  -     Pz     03P, 

(l-aiXl-aj)"^''     ''^■^ 
Sip  e  each  Sf  is  positive  definite  and  each  af  is  less  than  one,  each  of  the  terms  in 
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the  sum  is  negative,  unless /5_*  =  0  for  some  /. 

Let  us  now  consider  the  case  when  the  active  constraints  are  not  exactly 
satisfied,  i.e.  c*#0.  This  situation  does  not  arise  when  all  the  constraints  are  linear. 
However,  it  is  instructive  to  determine  how  far  the  results  for  the  linearly- 
constrained  case  can  be  extended  to  the  more  general  problem.  Although  dropping 
the  assumption  that  c''  =  0  means  that  the  theorem  no  longer  holds,  we  can  still 
show  that  while   y^v  m^y  be  an  ascent  direction,  ^Zfp^    will  still  be  a  descent 

direction. 

Theorem  4.7:  Suppose  that  A*  has  full  rank.  Then  S^fi'r  •  ^he  component  of 

i 

the  solution  of  the  PQP  that  is  orthogonal  to  the  range  space  of  A*  ,  is  a  descent 
direction  for  /  provided  it  is  not  zero.  In  referring  to  the  matrix-vector  product  //*x 
(4.4),  we  will  enclose  the  second  term  in  square  brackets  to  indicate  that  it  is 
irrelevant  in  the  event  that  the  diagonal  block  is  zero.  The  superscript  k  is  dropped 
as  above. 

Proof: 

gQP^    =     S     +     ffYlPy 

=  ^  +  {y^p^ 

^^Pz,  =  -Z\{g  +  [Kip,]) 

SQP,  =  8QP,  +   "I^^LP:, 

zIsqp^  =  (1  -  oii)Z2gQP^       as  above. 
=  (l-ai)Z^(^  +  [Kip,_]) 
=  (\-ai)Z2g  by  Lemma  4.5. 


Note  that  the  projected  gradients  are  the  same  as  for  the  case  when  Py  =  0.    The 
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ilue  of  Pv  fias  no  effect  on  the  computation  of  any  p,  because  of  the  orthogonality 
of  the  }',  and  Z,. 

Theorem  4.8:  Suppose  that  x*  is  a  feasible  point,  c''  =  0  and  A*  has  full  rank. 
Th^n  p*,  the  solution  of  the  MQP,  is  a  descent  direction  for  /,  provided  it  is  not 
zero . 

Proof:  As  defined  in  Section  (2.3).  /q  represents  the  working  set  at  the  start  of 
f.he  MQP  -^ubproblem.  /^  represents  the  working  set  at  the  end  of  the  MQP 
subPi-oblem,  ,4*  is  the  matrix  whose  rows  consist  of  the  transposed  gradients  of  the 
corresponding  constraints  and  t)''  are  the  QP  multipliers  corresponding  to  /j^.  We 
recaii  that  the  MQP  cannot  delete  any  of  the  constraints  that  were  intially  active. 
Therefore,  /f  =  ./q  U  Jn  where  7*;  are  the  inequality  constraints  that  are  active  at 
the  solution  to  the  MQP.  Thus,  at  the  solution  to  the  MQP,  we  have  that  either 
(df)jp  =  0  corresponding  to  one  of  the  initially  active  constraints,  or  ■r\'j  2  0  and 
(a|)jf  =  -c'^,  corresponding  to  a  constraint  that  was  added  to  the  working  set  by 
the  MQP.  Note  that  the  feasibility  of  v*"  implies  that  r*  >  0,  and  that  the  optimality 
cor.  itions  for  the  subproblem  imply  that  "n*  ^  0  for  all  j  i  J/^;. 

Let  r^  denote  the  vector  formed  by  choosing  those  components  of  c  which 
co(  espond  to  J^.    At  the  solution  to  the  MQP  we  have: 

g     +  H  p     -  AfT\ 

g'^p"  =  -p''h'p'  -  4\'  (4.7) 

^  0,  since  H''   is  positive  definite,   r*  =  0,  J  d  Jo, 
and  c),  T\j  ^  0,  jiJ^f,. 
Note  that  this  proof  does  not  apply  to  the  PQP,  since  some  ti*  may  be  negative. 

Theorem  4.9:  Suppose  that  v*  is  a  feasible  point,  t**  =  0  and  A*  has  full  rank 
Then  p*",  the  solution  of  the  IQP,  is  a  descent  direction  for  /.  provided  it  is  not  zero. 
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Proof:  It  follows  from  the  optimaiity  conditions  that  at  the  solution  to  the  IQP, 
■r\j  ^  0  for  all ;  €  Jf-.    Since  c*  2:  0,  the  proof  follows  as  for  Theorem  4.8, 

Global  Convergence  for  Linearly  Constrained  Problems 

Global  convergence  for  methods  which  modify  the  working  set  during  the 
major  iterations  is  strongly  dependent  on  the  stratgey  used  for  deleting  constraints 
from  the  working  set.  Various  strategies  for  modifying  the  working  set  which  do 
not  require  subspace  minimization  and  yet  avoid  zigzagging  have  been  studied. 
Byrd  and  Shultz  (1982),  Fletcher  (1981),  Kovacevic  (1975)  and  Ritter  (1980), 
among  others,  have  presented  such  methods  which  are  shown  to  be  globally 
convergent.  In  the  following,  we  shall  show  that  our  methods  fit  into  the 
framework  suggested  by  Byrd  and  Shultz.  Thus,  we  can  apply  their  results  to  prove 
that  our  EQP,  PQP  and  MQP  methods  are  globally  convergent. 

It  is  of  interest  to  note  that  McCormick  (1969)  proved  global  convergence  for 
an  algorithm  which  solves  the  bound-constrained  problem  by  "bending"  the  search 
direction  when  a  new  constraint  is  encountered.  In  fact,  our  PQP  algorithm  is  quite 
similar  to  this  method.  However,  McCormick 's  extension  of  his  bending  method  to 
the  general  linearly  constrained  problem  (1970)  differs  from  our  algorithm  since  his 
method  allows  for  moving  off  a  constraint  boundary  within  the  subproblem  if  the 
multiplier  estimate  is  negative. 

The  results  of  Byrd  and  Shultz  apply  to  algorithms  which  fit  into  the  following 
general  framework: 

0.  Given  x\  pick  A  ^  to  be  a  matrix  whose  rows  are  a 

linearly  independent  subset  of  the  constraints  satisfied  with 
equality  at  x^. 

Setk  =  1. 
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1.  Compute  Z^  such  that  4*Z*  =  0,  and  /?*  €  rangelZ''). 

2.  Compute  Lagrange  multiplier  estimates  X.*.    If  the 

dropping  rule  so  indicates,  then  drop  a  constraint  with 
a  negative  multiplier  estimate  and  update  Z"  and  p 
as  in  Step  1 . 

3.  Use  a  backtracking  line  search  along/;*  to  get  a  scalar  a*, 

and  augment  ,4*  if  a  constraint  is  hit  by  the  step  ci*p*. 

Set  -r*"^'  =  A*  +  a'^p''. 

4.  If  the  stopping  conditions  are  met,  then  stop 

5.  Set  /t  =  /t  4-  1  and  GOTO  Step  1. 

Our  EQP,  PQP  and  MQP  algorithms  follow  this  general  outline.  However, 
our  algorithms  update  A*^  within  the  QP  subproblem,  whereas  the  Byrd  and  Shultz 
framework  involves  modifying  the  working  set  within  the  line  search.  We  note  that 
when  the  full  step,  a*  =  1 ,  is  taken,  our  methods  reduce  to  the  Byrd  and  Shultz 
case.  When  the  full  step  is  not  taken,  a  simple  modification  to  our  algorithms  would 
ensure  that  they  satisfy  the  assumptions  made  by  Byrd  and  Shultz  If  we  do  not  take 
a  full  step  during  the  line  search,  we  can  immediately  modify  the  working  set  by 
removing  those  constraints  that  were  anticipated  to  be  added  by  the  QP  Allowing 
the  QP  to  add  constraints  to  the  working  set  and  modify  A",  followed  by  removal  of 
the  constraints  that  were  not  in  fact  encountered  by  the  line  search,  would 
essentially  be  the  same  as  adding  only  those  which  are  encountered  by  the  step 
a  p  . 

The  dropping  rule  proposed  by  Byrd  and  Shultz  is  simply  that  no  constraint  be 
dropped  from  the  working  set  if  the  previous  step  involved  adding  a  constraint, 
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unless  the  projected  gradient  is  zero.    This  dropping  rule  is  shown  to  guarantee 
global  convergence  under  the  following  assumptions: 

(1)  For  a  fixed  \x  >  0,  ^*^V*  ^  -tt||Z*'^*ll  lipll- 

(2)  For  any  subsequence  v  \  Z  ^  f  '  -  0  if  and  only  if  p  '  -  0. 

(3)  The  multiplier  estimates  \*  are  chosen  so  that  \f  x  '  ^  x  with  g  €/?(A  '  )  for  all 
;,  then  \*-X.*,  where  zero  multipliers  are  associated  with  the  inactive 
constraints. 

(4)  For  a  fixed  8  >  0,  whenever  a  constraint  is  dropped,  the  ratio  of  its  multiplier 
to  the  most  negative  multiplier  must  be  greater  than  8. 

(5)  The  steplength,  a*,  is  the  first  scalar  in  the  sequence  p^^(3^^  •  ■  •    satisfying 

fix")  -  fix^  +  ^Jp")  ^y^^g'^'p". 
where  3^  =  1  and  each  p^^^  ^  pp^  for  a  fixed  p  >  0  and  7^(0,1)- 

Note  that  the  least  squares  multipliers  satisfy  Assumption  3,  and  that  dropping 
the  constraint  with  the  most  negative  multiplier  estimate  is  consistent  with 
Assumption  4.  Assumption  5  clearly  can  be  satisfied  for  all  of  our  methods  as  well 
(see  Section  3.2),  since  all  of  our  methods  produce  descent  directions  for  /  .  It 
remains  to  be  shown  that  our  algorithms  satisfy  Assumptions  1  and  2. 

Global  convergence  of  the  EQP  algorithm  is  discussed  by  Byrd  and  Shuhz.  In 
the  following,  we  extend  their  results  to  our  PQP  and  MQP  methods. 

Lemma  4.10:  Let  p*  be  the  solution  to  the  PQP  at  .t^   Then 
|[p*||  <  K«||Z*^,,''||  for  some  k>0. 

Proof:  Consider  the  ith  step  of  the  PQP. 

^  iisni  wzfshp.w 

from  the  proof  of  Theorem  4.6  and 
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the  relationship  between  Z*  and  Zf 

Th'  s,  we  have 

\]^fp,;^\\^K\\Z\p,^'\\fo:a\\i. 
The    .emma  follows  immediately  from  the  fact  that  the  PQP  can  take  at  most  n 
step'- 

Lemma  4.11:  Assume  that  c*  =  0.  Let  p*  be  the  direction  produced  by  the 
PQP  subproblem  and  assume  that  af  >  0  for  some  /  Then  /?*  satisfies  Assumption 
1. 

Proof: 

-g'^V*  =  -iQpp''  from  ;4.5) 

=  +  i2\hy  fa\B\-\z'(g'Qp)  +  iZ'jg'Qpya\Br{ZUhp)  +    '  '  ' 
=  R'[{Zfg'Q,fB\-\.zWbp)  +  (l-"f)a2*(Zf^e/'/«2"'(Zf4p) 

+  (l-Sf)(l-a2*)53*(Z3*'4^/sr'(Zf4^_)  +    •  •  • 
by  Theorem  4.6. 

>  a^WBVr'  ||zf4^JP  +  (l-5f)S*||S^||-^  \\Z'2g'Q,f  +    ■  •  • 

ilstllllfir'l! 

"  1^'''^*''^^^'"  '^*''  ^^  Lemma  4.10,  for  ?  >  0. 
Note  that  if  af  =  0,  it  can  be  replaced  by  the  first  nonzero  af. 

,\lthough  for  Lemma  4.11  we  have  assumed  that  a,'"  #  0  for  some  ;,  this 
assumption  is  not  a  restriction.  If  af  =  0  for  all  ;.  then  the  subproblem  has  simply 
picked  up  many  constraints  without  taking  a  step.  Since  constraints  have  been 
added  to  the  working  set,  the  dropping  criterion  will  net  be  met  and  the  next 
iteration  will  start  with  the  projected  gradient  equal  to  zero. 
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The  MQP  search  direction  is  produced  by  taking  a  PQP  step  and  then  possibly 
moving  off  one  of  the  constraints  that  the  PQP  added  and  continuing  in  the  new 
space.  As  discussed  in  Section  2.4,  the  MQP  moves  off  a  constraint  added  by  the 
PQP  only  when  there  is  another  constraint  that  is  more  binding  at  the  solution  to  the 
subproblem  (see  Figure  2.1).  Therefore,  the  MQP  will  never  retrace  the  steps  of 
the  PQP  in  the  opposite  direction.  Assumption  1  thus  holds  for  our  MQP  method  as 
well. 

To  see  that  Assumption  2  holds  for  both  our  PQP  and  MQP  methods,  we 
consider  the  geometry  of  the  subproblem  (see  Figure  2.1).  It  is  immediately 
obvious  from  the  definition  of  p*  for  the  PQP,  that  if  Z*"^^*  -  0,  p*  -  0  as  well.  For 
the  MQP,  p*  starts  out  as  the  PQP  step,  but  then  may  include  some  extra  steps  by 
moving  off  some  constraints  that  had  just  been  added.  It  is  obvious  from  the 
geometry  that  moving  off  any  constraint  that  has  been  added  to  the  working  set 
cannot  cause  the  MQP  direction  to  move  "backward"  along  the  PQP  direction. 
Therefore,  it  follows  for  the  MQP  direction  as  well,  that  if  Z*'^*  -  0,  p*  -  0. 

For  the  reverse  implication,  we  first  assume  that  af  #  0.  Again,  from  the 
geometry  of  the  subproblem,  we  can  see  that  both  the  PQP  and  MQP  steps 
pf,  I  >  1,  cannot  point  back  in  the  direction  -p*.  Thus,  if  p*  -0,  it  follows  that 
ajZ^i-O  and  from  the  definition  of  pi  for  either  method,  Z*'^*  -  0.  If  we 
remove  the  assumption  that  a*  +  0,  it  may  indeed  be  the  case  that  p*  is  close  to 
zero  but  Z*  ^*  was  not  near  zero.  However,  since  this  means  that  some  constraint 
was  added  within  the  QP,  the  dropping  condition  will  not  be  met,  and  the  projected 
gradient  at  the  start  of  the  next  QP  will  be  near  zero. 

We  can  now  state  the  global  convergence  result: 

Theorem  4.12:  Suppose  /iR"  -  R,/  €  C^R"),  and  x'  is  a  cluster  point  of  the 
iterates  generated  by  either  our  EQP,  PQP  or  MQP  algorithms  using  a  dropping 
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rule  dtisfying  the  Byrd  and  Shultz  full  step  dropping  condition,  and  either  ot*^  =  1 
for  all  k,  or  we  use  the  simple  modification  mentioned  above.  Then  if  the  active  set 
at  x'  ::i  linearly  independent,  x    is  a  Kuhn-Tucker  point. 

Proof:  The  proof  follows  immediately  from  Byrd  and  Shultz's  Theorem  1  and 
the  fact  that  our  methods  satisfy  all  of  the  Byrd  and  Shultz  assumptions. 

The  next  theorem,  also  due  to  Byrd  and  Shultz,  establishes  that  if  the  solution 
satisfies  strict  complementarity  and  the  sufficient  conditions  for  a  strict  local 
minimum,  then  the  working  set  settles  down  and  does  not  have  to  be  repeatedly 
modified. 

Theorem  4.13:  Suppose  that  the  point  i-  is  a  cluster  point  of  a  sequence 
generated  by  our  EQP,  PQP  or  MQP  methods,  under  the  assumptions  made  for 
Theorem  4,12.  that  c'  is  a  Kuhn-Tucker  point  such  that  gix")  +  A*  \'  =0  with  A* 
of  full  rank,  X*  >  0,  ali  other  constraints  are  strictly  satisfied  and 
Z(x')^W(x\\')Z(x^)  positive  definite.  Then  the  sequence  converges  to  .v*  and  the 
active  set  remains  constant  for  all  k  sufficiently  large 

Proof:  The  proof  follows  from  Byrd  and  Shuhz's  Theorem  2. 


4A  A  PQP-type  Method 

iince  we  are  able  to  prove  that  the  direction  obtained  by  solving  the  PQP  is  a 
descent  direction  in  the  case  that  the  active  constraints  are  exactly  satisfied 
(Theorem  4,6),  it  is  instructive  to  consider  a  PQP-type  method  for  which  this 
condition  would  always  bold.  The  proposed  method,  which  we  call  PQP2,  would 
start  each  QP  subproblem  with  no  acti\  e  constraints.  The  QP  would  then  add 
constraints  to  the  working  set  as  they  are  encountered.  Since  this  method  would 
start  each  QP  with  a  nul.  'vorking  set,  the  problem  of  which  constraints  to  delete  is 
solved  by  simply  dropping  them  all.    In  the  neighborhood  of  a  solution,  it  is  to  be 
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expected  that  only  the  constraints  in  7*  will  be  picked  up. 

Clearly,  this  method  is  highly  inefficient.  We  expect  the  working  set  at  one 
iteration  to  be  fairly  close  to  the  working  set  of  the  previous  iteration.  It  is 
advantageous,  therefore,  for  the  PQP  method  to  start  the  QP  subproblem  with  the 
working  set  of  the  previous  iteration.  Obviously,  by  disregarding  the  working  set  of 
the  previous  iteration,  the  proposed  method  will  entail  a  considerable  amount  of 
overhead  in  picking  up  the  active  constraints  at  the  start  of  each  subproblem.  More 
importantly,  by  starting  with  a  null  working  set,  this  method  essentially  resets  the 
projected  Hessian  approximation  at  each  iteration.  Nevertheless,  it  is  instructive  to 
examine  the  behavior  of  such  a  method,  in  order  to  properly  understand  the  role  of 
the  active  set  selection. 

We  are  able  to  show  that  when  the  method  reaches  a  solution,  assuming  that 
the  gradients  of  the  active  constraints  are  linearly  independent,  the  algorithm  will 
indeed  identify  the  true  active  set.  However,  we  demonstrate  that  this  method  can 
experience  difficulty  in  moving  away  from  a  non-optimal  vertex.  The  difficulty  can 
be  resolved  only  by  considering  the  constraint  multipliers.  Our  results  emphasize 
the  crucial  role  of  correct  active  set  determination  and  the  necessity  for  accurate 
multiplier  estimates. 

At  the  start  of  each  iteration  of  PQP2  there  are  no  constraints  in  the  working 
set.  Since  the  gradient  can  be  expressed  as  a  linear  combination  of  some  subset  of 
the  constraint  gradients,  those  constraints  should  be  added  to  the  working  set.  For 
simplicity,  assume  that  the  problem  is  not  degenerate,  so  that  the  active  set  is  well- 
defined. 

The  initial  search  direction  is  given  by  p  =  -S*  ^*.  Clearly,  g''  p  <  0. 
Assume  that  there  is  some  constraint,  Cj,  such  that  Qj  p  <  0.  A  Taylor  Series 
expansion  of  the  ;th  constraint  is  given  by 


It  is  clear  that  ?f  cfc.r*)  was  exactly  satisfied,  i.e.,  Cjix'')  =  0.  a  step  along  p  would 
cause  it  to  be  violated.  In  fact,  the  algorithm  determines  the  steplength  by 
examining  the  sign  of  afp.  Hence  the  method  will  immediately  pick  up  some  such 
constraint  for  which  a)  p  <  0.  '    ■ 

As  constraints  are  added  to  the  working  set,  the  dimension  of  Z'^  is  decreased. 
The  search  direction  sben  become^;  -Zfsf'^zfg.  Since  Zf  by  construction  is 
orthogonal  to  the  active  constraint  gradients,  a'j  p  <  0  only  for  some  constraint  Cj 
that  is  not  yet  in  the  working  set.  This  new  constraint  will  then  be  added  to  the 
working  set.  Thus,  at  the  solution  the  constraints  which  define  the  gradient  will  be 
included  in  the  working  set. 

The  Allowing  example  demonstrates  the  failure  of  the  PQP2  method. 
Consider  the  problem 


min   xi 


s.t.    —xT  T  X2  —  1^0 

x}  +  ji'2  -  2  >  0 
-  .r,  +  I  >  0. 

The  last  consiraint  >s  added  to  limit  the  feasible  region  ao  that  the  solution  is 
not  unbounded.  The  first  two  constraint  boundaries  are  parabolas  which  intersect  at 
the  points  (v2,„0')^  and  --^2,0)'^.  The  second  parabola  has  a  higher  apex  and  a 
steeper  slope  see  Figure  4.1).  Therefore,  when  the  constraints  are  approached 
from  the  feasible  side  in  either  of  the  lower  quadrants  c;  is  always  encountered 
first  The  algorithm  will  proceed  by  stepping  to  the  first  constraint  and  moving 
along  that  linearized  constraint  until  it  reaches  the  second  one.  Moving  m  this 
manner,  the  i^^iates  will  eventually  converge  to  the  point  of  intersection  of  the  two 
constraints.    At  this  vertex  both  constraints  are   active      However,   the  objective 
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function  can  still  be  decreased  by  moving  off  the  first  constraint  while  remaining  on 
the  second. 

At  this  stage  the  method  can  be  said  to  have  failed,  since  it  has  converged  to  a 
non-optimal  point.  Moreover,  we  will  demonstrate  that  this  algorithm  will  not  be 
able  to  recover  by  moving  away  from  that  point.  Any  method  which  is  to  succeed 
when  applied  to  a  problem  of  this  nature  must  include  a  mechanism  for  moving 
away  from  non-optimal  vertices.  Such  a  move  must  be  based  upon  a  choice  of 
which  constraints  to  keep  active  and  which  to  allow  to  become  inactive.  The  PQP2 
algorithm  does  not  incorporate  any  active  set  strategy.  This  method  simply 
considers  none  of  the  constraints  active  until  they  are  encountered,  but  once  they 
are  reached  they  remain  active  for  the  duration  of  the  subproblem.  For  the  problem 
under  consideration,  the  behavior  of  the  algorithm  hinges  upon  which  of  the  two 
constraints  at  that  vertex  is  picked  up  first.  If  the  first  constraint  is  made  active, 
then  the  second  is  encountered  immediately  and  no  progress  can  be  made.  If, 
however,  the  second  constraint  is  picked  up  first,  then  the  method  can  step  along  the 
search  direction  in  the  null  space  of  the  second  constraint  gradient  without  reaching 
the  first  constraint. 

We  now  describe,  in  detail,  the  performance  of  the  algorithm  on  this  problem. 
The  gradient,  g,  is  (1,  0)^  and  the  matrix  of  transposed  constraint  gradients.  A,  is 
given  by 

A  = 

The  solution  is  at  x'  =  (  —  2,  - 1)^. 

Let  x^  =  (3,  -1)^.  The  PQP2  would  proceed  as  follows:  The  initial  search 
direction,  p,  is  along  the  negative  gradient,  (-1,0)^.  As  we  search  along  p,  the 
linearization  of  the  constraints  provides  the  approximation 
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cix^  +  p)  =  c(x^)  +  A°p. 
The'cfore,  the  predicted  value  of  the  constraints  is  given  by 


:(.v"  +  p)~ 


'2.5' 

3      1  ' 

T 

6.0 

+ 

6     1 

2.5. 

0    -1 

r-1 

0 


-.5' 
0      . 
2.5, 
Thus,   .1  unit   step  along  p    would   cause  the   first  constraint  to  be  violated      The 

algorithm  therefore  steps  along  p  until  it  reaches  the  linearized  constraint  boundary 

for  fi-    The  next  search  direction  is  in  the  tangent  space  to  the  gradient  of  c,,  which 

is  given  by  (-.0167,  .0500)^^.    Stepping  along  this  direction,  the  second  constraint  is 

encountered,  and  the  PQP  is  terminated  giving  t'  =  (2.150,  -.9500)^. 

The  nonlinear  functions  are  reevaluated  outside  the  PQP  and  the  constraints 
are  reapproximated  by  their  linearizations  at  the  new  point.  Once  again,  a  step 
along  the  negative  gradient  causes  c-  to  be  encountered  first,  and  then  a  step  can  be 
taken  by  remaining  on  that  linearized  constraint.  The  algorithm  steps  to 
x^  -  (1.834,  -.6319)^  A  similar  process  is  repeated  at  each  iteration,  producing 
the  iterates  v^  =  (1.584.  -.2230)^.  v^  =  (1.4233,  0)^.  v^  =  (1.414,  0)^  which 
cor  verge  to  .i  =  i^l,  0)   . 

Both  constraints  cj  and  ci  are  active  at  v.  However,  v  is  not  an  optimal  point, 
since  the  solution  to 

A^K  =  g 
is  given  by  \i  =  ( 1— ,  ~;~")   >  v^hich   indicates  that  the   first  constraint  can  be 

dropped  This  method,  however,  does  not  allow  for  dropping  a  constraint  within 
the  PQP.  The  algorithm  would  simply  delete  all  the  constraints  from  the  working 
set  and  reenter  the  PQP. 

The  algorithm  will  fail  at  this  stage  if  Cj  is  again  picked  up  first     Since  any  step 
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along  ci  will  violate  C2,  the  second  constraint  will  have  to  made  active  immediately 
and  no  progress  will  be  made  towards  the  solution.  However,  were  C2  to  be  added 
to  the  active  set  first,  a  step  could  be  taken  without  encountering  any  other 
constraints,  producing  the  iterate  (1.3031,  .3143)^. 

The  difficulty  this  method  experiences  in  moving  away  from  a  non-optimal 
vertex  can  be  resolved  only  by  considering  the  Lagrange  multiplier  estimates. 
Examples  of  this  nature  are  characterized  by  having  a  choice  of  constraints  to  be 
made  active.  Rather  than  arbitrarily  choosing  any  one,  the  selection  of  constraint  to 
be  picked  up  could  be  based  upon  the  multiplier  estimates.  Of  course,  this  raises 
the  question  of  whether  such  estimates  are  reliable.  In  a  case  such  as  the  one 
described  above,  where  the  PQP  is  started  at  a  vertex,  (or,  indeed,  any  nonoptimal 
stationary  point)  the  muhiplier  estimates  would  be  based  on  accurate  constraint 
values,  since  the  constraints  and  their  gradients  are  evaluated  immediatley  prior  to 
the  start  of  the  PQP.  In  any  other  case  in  which  the  PQP  steps  to  a  vertex,  we  can 
allow  the  algorithm  to  proceed  as  usual.  Should  the  PQP  then  terminate  as  a  result 
of  having  too  many  active  constraints,  at  the  next  iteration,  having  reevaluated  the 
constraints,  either  some  will  no  longer  be  active,  or  the  PQP  will  be  started  at  a 
vertex  and  the  multiplier  information  can  be  used. 

Such  a  method  would  indeed  guarantee  that  eventually  the  algorithm  would 
move  off  a  constraint  with  a  negative  multiplier.  However,  even  in  an  example  such 
as  the  one  considered  above,  the  algorithm  takes  many  short  steps  until  it  reaches 
the  stage  at  which  the  QP  is  started  at  a  vertex.  Therefore,  although  a  method  of 
this  nature  would  have  interesting  theoretical  properties,  in  practice  it  would  not  be 
competitive. 
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Chapter  5 
Numerical  Results 


5.1  Test  Problems  and  Results 

The  results  were  obtained  by  running  the  NPSOL  package  Version  2.0,  and  a 
Fortran  implementation  of  our  methods,  on  a  VAX-11/780  at  the  Courant 
Mathematics  and  Computing  Laboratory.  Double  precision  arithmetic  was  used 
throughout. 

As  discussed  in  Section  3.2,  there  is  no  merit  function  with  which  our  methods 
will  perform  uniformly  well.  We  begin  by  presenting  two  problems  which 
demonstrate  how  either  the  exact  penalty  merit  function  or  the  augmented 
Lagrangian  merit  function  used  in  NPSOL  can  degrade  the  performance  of  our 
methods,  while  using  the  other  merit  function  on  the  same  problem  would  not  cause 
any  difficulty. 

The  first  problem,  presented  in  Chamberlain  et.  al.  (1982),  demonstrates  the 
"Maratos  effect"  which  can  arise  from  the  use  of  the  exact  penalty  merit  function. 
The  second  problem  shows  how  use  of  the  NPSOL  merit  function  may  cause  our 
methods  to  take  many  extra  steps.  This  problem  also  demonstrates  the  sensitivity  of 
the  exact  merit  function  to  the  correct  choice  of  the  penalty  parameter,  p.  The 
results  reported  are  for  p  =  10.  In  fact,  since  X*  =  (-1,  -1)^,  the  same  results  can 
be  obtained  for  any  p  ^  1.  However,  setting  p  =  .9  causes  all  of  our  methods  to 
fail,  since  the  search  directions  are  no  longer  descent  directions  for  the  exact  penalty 
merit  function. 
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The  algorithms  run  are  described  in  detail  in  Chapter  2.  For  each  run,  we 
indicate  which  merit  function  was  used,  and  for  the  exact  penalty  merit  function  the 
value  of  the  penalty  parameter  p.  Unless  otherwise  indicated,  the  Byrd  and  Shultz 
dropping  rule  was  not  enforced,  the  tolerance  level  for  dropping  constraints  was  set 
to  1.^  —  3  and  definitions  (a)  and  (c)  of  Section  3  4  were  used  for  s'^  and  v*  in  the 
qua.M-N'ewton  BFGS  update.  The  matrix  D,  used  in  forming  the  approximation  to 
the  full  Hessian  matrix  from  the  projected  matrix,  as  in  (2.12),  was  set  to  the 
identit\  matrix. 

The  convergence  criteria  are  defined  as  in  the  \PSOL  package.  There  are 
essennally  four  conditions  tested: 

1        ilzVll^fo/i 

2.  NVIOL  =  0 

3.  X^  >  toll 

4        a*iip*||  <  toll 
where 

toll  =  (10.0€„,;,)''-max(l  +  [f(.v*)|,:|^>/j||), 
toll  =  i\O.Oi„,,)''H\  +  \[x'\\) 
and   ./^^  refers  to  the  gradient  with  respect  to  the  free  variables.    Convergence  is 

de;ii  ed  as  either  the  satisfaction  of  all  four  conditions,  or  of  conditions  2  and  3 

alo".g  with 

5.       l|zVll^emc.max(l  +  l/(.r*)|,||^*;j||). 

The  following  quantities  are  reported  for  each  method: 

NMJ  -   the  number  of  "major"  iterations; 

NMN  -  the  number  of  "minor"  iterations; 

NF   -  the  number  of  function  calls  (i.e.,  the  number  of  evaluations  of  f,  g,c  and 
A). 
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An  entry  of  stars  indicates  failure;  a  dashed  line  indicates  that  the  corresponding 
combination  of  method  and  parameter  does  not  apply.  For  each  problem  we  list 
below  the  number  of  variables,  n,  the  number  of  bound  constraints,  mb,  the  number 
of  linear  constraints,  ml,  and  the  number  of  nonlinear  constraints,  mn. 

The  number  of  function  calls  is  determined  by  the  number  of  trial  points  the 
merit  function  must  consider  before  the  new  iterate  is  accepted.  At  best,  the  merit 
function  always  accepts  the  unit  step,  and  NF  =  MMJ+l.  However,  for  the  EQP, 
and  sometimes  the  PQP  and  MQP  methods,  the  algorithm  may  take  a  "zero"  step, 
i.e.,  add  some  constraints  to  the  working  set  without  making  a  move.  In  such  a 
case,  the  number  of  major  iterations  may  be  larger  than  the  number  of  function 
calls. 

Problem  5.1.  Chamberlain  et.  al.  (1982).  n  =  2,  mb  =  0.  ml  =  0.  mn  =1. 
The  matrix  Z  has  1  column  at  the  solution. 

/=  -xi  +  \0{xl  ^  x\-\) 

s.t.   C]  =  .ti  +  jrl  -  1  =  0 
.v°  =  (.8.  .6)^ 

/*  =  -1 

x'  =  (1.0,0.0)^ 

Problem  5.2.  Mukai-Polak  (1978).  n  =  6.  mb  =  1 .  ml  =  0,  mn  =  2.  The 
matrix  Z  has  3  columns  at  the  solution. 

f=(Xi-  X^f   +    (X2  -   X5)2   +    (.V3   -   X6)2 

s.t.    ci  =  xj  +  .r^  +  ;r|  <  5 
C2  =  xl  +  (Xi  -  3)2  <  1 
4  <  .V6  <  8 
x°  =  (1.0,  1.0,  1.0,  3.0,0.0,  0.5)^ 
/*  =  5 
X*  =  (1.0,0.0,  2.0,  2.0.0.0,4.0)^ 
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Table  5.1 


Chamberlain 

Pe 

Method 

Ps' 

p  =  10 

p  =  0.9 

p  =  0.1 

NPSOL 
IQP 

6,7,7 
5,5,6 

31,31,146 

18,19,68 

17,18,63 

MQP 

5,5,6 

31,31,146 

18,19,68 

17,18,63 

PQP 

5,5,6 

31,31,146 

18.19,68 

17,18,63 

EQP 

5.5.6 

31.31,146 

18,19,68 

17,18,63 

Table  5.2 


Mukai-Polak 

Method 

^.v 

Pe 

p  =  10 

p  =  0.9 

p  =  0.1 

NPSOL 

IQP 

MQP 

PQP 

EQP 

9,12,15 
27,30,53 
26,29,51 
26,29,51 
10,11,14 

11,14,14 
10,13,13 
10,13,13 
11,12,13 



Incorrect  choice  of  merit  function  can  cause  our  methods  to  fail  or  to  reject  the 
full  step  a  =  1   even   near  the  solution      Rejection  of  the  full  step  has  a  doubly 
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negative  effect  on  the  performance  of  these  methods.  Firstly,  if  the  full  step  is  not 
taken,  more  of  these  shorter  steps  must  be  taken,  causing  the  number  of  major 
iterations  to  increase.  Secondly,  the  NPSOL  package,  and  our  modifications  of  it, 
use  the  working  set  from  the  end  of  the  previous  subproblem  as  the  intial  working 
set  for  the  next  iteration.  Therefore,  if  the  full  step  is  not  taken,  it  may  be  the  case 
that  constraints  that  are  included  in  the  working  set  because  they  were  added  by  the 
previous  QP  are,  in  fact,  not  on  their  bounds  because  the  full  QP  step  was  not 
taken.  In  that  case,  at  the  next  iteration  the  feasibility  phase  may  cause  the 
algorithm  to  take  extra  steps  in  order  to  move  onto  those  constraints,  thereby 
increasing  the  number  of  minor  iterations. 

For  these  reasons,  we  are  interested  in  studying  the  performance  of  our 
methods  without  being  concerned  with  the  possible  adverse  affects  of  a  bad  choice 
of  merit  function.  The  minimum  ly,  and  /j  problems,  discussed  in  Section  4.2, 
provide  us  with  a  class  of  problems  for  which  a  "natural"  merit  function  is  available. 
Note  that  although  our  implementations  are  a  reasonable  way  to  solve  minimum  l,c 
problems,  they  give  a  very  inefficient  way  of  solving  minimum  /j  problems  because 
of  the  extra  variables  introduced.  See  Overton  (1982)  for  a  discussion  of  efficient 
implementations  of  minimum  /^  methods. 

The  test  problems  listed  below  can  be  found  in  Gill  and  Murray  ( 1976).  Many 
of  these  problems  are  also  discussed  by  Watson  (1979)  and  More  et.  al.  (1981). 
These  problems  were  originally  used  for  testing  minimum  least  squares  methods, 
but  they  are  of  interest  for  testing  other  data-fitting  methods  as  well.  Note, 
however,  that  the  minimizer  and  minimum  values  reported  in  the  literature  for  the 
least  squares  case  usually  differ  from  the  optimal  values  for  the  /j  and  /^  norms. 

For  each  of  the  test  problems,  we  list  below  the  number  of  variables,  n,  the 
number  of  functions,  m,  the  initial  point,  the  minimizer  and  minimum  objective 
value  obtained  by  NPSOL  for  both  the  minimum  /j  and  minimum  /x  norm.    For  the 
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problems  on  which  NPSOL  failed,  the  resuhs  listed  are  those  obtained  by  our  IQP 
methcd.  Our  methods  found  a  different  solution  from  that  found  by  NPSOL  for 
both  the  minimum  U  and  the  minimum  l^  norms  on  Problem  Chebyquad.  The 
solutions  found  by  our  methods  are  listed  following  the  NPSOL  results.  Note  that 
for  minimum  /j  problems,  the  equivalent  nonlinear  programming  problem  has 
n  =  n  ^  in  variables  and  mn  =  2in  constraints  Similary,  the  nonlinear 
programming  problem  for  the  minimum  I „  problem  has  «  =  n  -f  1  variables  and 
mn  =  2w  constraints. 

Problem  5.3.    Woods,    n  ■■=  4.     w  =  7. 

r'^  =  (2.0,  2.0,  2.0,  2.0)^ 

l^:  x'  =  (1.0000,  1.0000.  1.0000,  1  0000)'' 

/*  =  .89099^-07 
/j:   x'  =  (1.0000,  1.00(X),  1.0000,  1.0000)^ 

/*  =  0.0000 
Problem  5.4.    Singular,    n  =  4,     m  =  4. 

.x°  =  (2.0,  2.0,  2.0,  2.0)^ 

/^:   X*  =  (0.60726e-08,  -  0.60726^-08.  0.654;4£'-09, 
0.65414^-09)'" 
/*  =  0.0000 
/i:   .r'  =  (0.29440^-03,0.29440^-04,0.15152^-02, 
0.15152^-02)'" 
/•  =  0.0000 
Problem   5.5.    Rosenbrock.   n  =  2,     m  =  2. 

x°  =  (2.0,  2.0)^" 

/^:   x'  =  ri.OOOO,  1.0000)'" 

/*  =  0.00000 
/i:   X*  =  (1.0000,  1.0000)'" 

/  =  0.0000 
Problem  5.6.    Chebyquad.    n  =  8.     m  =  8. 

.v°  =  ( \/n+  1 .  2AT+  1 ,       ■  ■  ,  n/n+  1 )'" 
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U:   x'  =  (0.62000^-01,0.25867,0.29097,0.511897,0.511897, 

0.73992,  0.812691,  0.961794)^ 
/*  =  0.242579^-01 

x'  =  (0.686159^-01,  0.21318,  0.29142,  0.50000,  0.50000, 
0.70858,0.78682,0.93138)^ 

/*  =  0.42939^-01 
/j:   X*  =  (0.1001,  0.25125,  0.317180.5000,  0.5000,  0.68282, 

0.74875,  0.89902)'" 

/*  =  0.15764 

x'  =  (0.620246^-01,  0.21519,  0.25425,  0.50000,  0.50000, 
0.74575,0.78481,0.93797)^ 

/*  =  0.952535e-01 

Problem  5.7.    Helix,    n  =  3,     w  =  3. 

.r°  =  (2.0,  2.0,  2.0)^ 

U:   x'  =  (1.0000,  -0.32378^-10,  0.0000)^ 

/*  =  0.0000 
/^.   x'  =  (1.0000,  -0.30811^-09,0.0000)'" 

/*  =  0.0000 
Problem  5.8.    Large  Residual  (Davidon  2).   «  =  4,     m  =  20. 

x°  =  (25.0,5.0,  -5.0,  -1.0)'" 

l^.   /  =  (-12.2437,  14.0218,  -0.45151,  -0.10517^-01)'" 

/  =  115.7064 
/^.   ,/  =  (-10.2236,  11.9084,  -0.45804,0.58032)'" 
/*  =  903.2343 
Problem  5.9.    Cragg  and  Levy.    «  =  4,     m  =  5. 

x^  =  (2.0,  2.0,  2.0,  2.0)'" 

U:   .X*  =  (0.117127^-04,  1.0000,  1.0000,  1.0000)'" 

/•  =  0.0000 
l^.    /  =  (-0.2762^-04,  1.0000,  1.0000,  1.0000)'" 

/•  =  0.0000 
Problem  5.10.    Beale.    «  =  2,     m  =  3. 

,v°  =  (2.0,2.0)'" 
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U:   x'  =  (3.0000,  5.0000)^ 

/'  =  0.0000 

/j:   x'  =  (3.0000,  5.0000)'" 

/'  =  0.ll3956e-07 

Problem  5.11.   Branin.    «  =  2,     m  =  2. 

r°  =  (2.0,  2.0)'' 

/„:   x"  =  (0.0000,  O.OOOO)'" 

/'  =  0.0000 
ll.   x'  =  (-0.1916«?-08,  0.1916t--08)^ 

/*  =  0.0000 
Problem  5.12.    Freudenstein  and  Roth,    n  =  2.    m  =  2. 

x^  =  (15.0,  -2.0)^ 

l^:  .x'  =  (11.41278,  -0.89681)'" 

/*  =  4.948952 
/;:   x'  =  (15.51727,  -0.8968)'" 

/*  =  9.89790 
Problem  5.13.   Bard,   n  =  3,     m  =  \5. 

x^  =  (1.0,  1.0,  1.0)'" 

/»:   x'  =  (0.53469^-01,  1.53425,  1.96574)'' 

/*  =  0.5081633^-01 
/i;   x'  =  (0.10094,  1.5251,  1.97211)'" 

/*  =  0.12434 
Problem  5.14.    Kowalik  and  Osborne,    n  =  4,    m  =  ll. 

x°  =  (0.25,  0.39,  0.415,  0.39)'" 

U:   X*  =  (0.18463,  0.101515,  0.11935^-01,  0.11177)'" 

/"  =  0.80876^-02 
/j:   x'  -  (0.19337,0.19376,0.10893,0.13973)^ 

/•  =  0.38768e-01 
Problem  5  15.   Osborne  \.    n  =  5,     m  =  33. 

x°  =  (0.5,  1.5,  -1.0,  .01,  .02)'" 

U:   x'  =  (0.36636,  1.4052,  -0.92482,  0.11309^-01,  0.25575^-01)'" 
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/'  =  0.278179^-02 
/j:   x'  =  (0.37706,  2.1925,  -1.7255,0.13325^-01,  0,21288^-01)^ 
/*  =  0.293912^-01 

Problem  5.16.    Engvall  1.^=5,     m  =  5. 

a:°  =  (2.0,  2.0,  2.0,  2.0,  2.0)^" 

U:   x'  =  (0.0000,  0.0000,  1.0000,  3.4747,  5.1465)^ 

/•  =  0.0000 
Zj:   all  methods  failed 
Problem  5 .17.   Zangwill  2.    «  =  3,     m  =  3. 

x°  =  (2.0,  2.0,  2.0)^ 

U:   X*  =  (0.0000,  0.0000,  O.OOOO)'" 

/*  =  0.0000 
/j:   x'  =  (0.0000,  0.0000,  0.0000)^ 

/*  =  0.0000 

Problem  5.18.    Madsen.    n  =  2,     m  =  3. 

x°  =  (0.6,  -0.8)^ 

/„:   x'  =  (0.45330,  -0.90659)^ 
/*  =  0.61643 

/i:  jc'  =  (0.0000,  -0.181178^-02)^ 
/*  =  1.0000 
Tables  5.3  -  5.5  compare  the  performance  of  our  methods,  using  the  objective 
function  itself  as  the  merit  function,  with  NPSOL,  which  uses  an  augmented 
Lagrangian  merit  function.  For  the  minimax  problem,  we  ran  both  versions  of  the 
subproblem  discussed  in  Section  4.2.  As  before,  the  three  values  indicated  for  each 
method  represent  the  number  of  major  iterations,  the  number  of  minor  iterations 
and  the  number  of  function  calls.  An  entry  of  stars  in  the  table  indicates  failure. 
The  Byrd  and  Shultz  dropping  rule  was  not  enforced  and  the  tolerance  level  for 
deleting  constraints  was  set  to  l.e-3.  In  fact,  all  of  the  problems  were  also  run 
using  the  default  tolerance  level  (see  Section  3.5),  and  produced  nearly  identical 
results.    We  found  that  using  definition  (c)  of  Section  3.4  for  >•*  in  the  quasi-Newton 


-96- 

update  did  not  perform  well  for  the  data-fitting  problems,  possibly  because  the 
gradient  of  the  objective  function  remains  constant.  Instead,  we  use  definition  (0 
for  these  problems. 

As  evidenced  by  the  results  shown,  our  methods  generally  converge  somewhat 
faster  than  NPSOL  on  this  class  of  problems.  We  note  that  for  these  problems  we 
use  a  different  merit  function  than  the  one  used  in  NPSOl  ,  However,  for  problems 
on  which  NF  is  about  the  same  as  \MI .  the  merit  function  did  accept  steplengths  of 
one,  and  the  different  results  can  be  attributed  to  the  differencev  in  the  second 
derivative  information.  A  comparison  of  Tables  5.3  and  5.4  indicates  that  the  two 
formulations  of  the  subproblem  for  the  minimax  problems  can  produce  different 
results  but  that  neither  one  is  uniformly  better  than  the  other. 

The  column  of  the  tables  headed  "nz"  indicates  the  number  of  columns  of  Z  at 
the  solution,  which  is  the  dimension  of  the  projected  Hessian  matrix.  Many  of  these 
problems  are  characterized  by  having  a  large  number  of  constraints  active  at  the 
solution,  and  hence  the  reduced  space  is  of  much  lower  dimension  We  note  that 
although  a  zero  value  of  nz  indicates  that  the  projected  Hessian  is  null  at  the 
solution,  it  generally  is  the  case  that  there  is  some  matrix  of  small  dimension 
recurred  during  intermediate  stages  of  the  computation,  until  the  working  set  is 
correctly  identified. 

As  discussed  in  Section  1.4,  methods  which  update  an  approximation  to  the  full 
second  derivative  matrix  often  introduce  an  augmenting  term  in  order  to  force  the 
approximating  matrix  to  remain  positive  definite.  These  problems  demonstrate  the 
ill-conditioning  of  the  approximating  matrices  which  can  result  from  adding  this 
augmenting  term.  The  estimated  condition  number  of  the  NPSOL  Hessian 
approximation  for  the  Woods  minimax  problem  was  on  the  order  of  10  ,  for  the 
Helix  minimax  problem.  10^  and  for  the  Cragg  and  Levy  minimax  problem,  10  , 
However,  for  none  of  the  minimax  problems  did  the  estimated  condition  number  of 
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the  projected  Hessian  approximation  recurred  by  our  methods  grow  larger  than  10^. 
We  note  that  for  some  problems,  our  methods  indicate  a  very  large  number  of 
function  calls,  although  the  number  of  major  and  minor  iterations  may  be  small.  In 
most  cases,  this  phenomenon  occurred  on  problems  for  which  some  Lagrange 
multipliers  at  the  solution  were  near  zero.  In  such  cases,  the  methods  may  have 
essentially  converged  earlier.  At  the  last  step,  the  merit  function  may  consider 
many  trial  points  in  an  attempt  to  find  the  next  iterate,  causing  the  number  of 
function  calls  to  increase. 
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Table  5.3 


Minimum  /^  problems  -  Version  1 

Problem 

nz 

NPSOL 

IQP 

MQP 

PQP 

EQP 

Woods 

0 

32,88,83 

6,16,7 

7,12,7 

7,12,7 

10,11,9 

Singular 

1 

31,51,33 

33,10,160 

33,61,347 

33,61,347 

15,30,48 

Rosenbrock 

1 

24,29,38 

8,22,11 

9,11,18 

9,11,18 

11,11,22 

Chebyquad 

5 

28,57,33 

6,11,8 

6,11,8 

6,11,8 

10.11,15 

Helix 
Large  Res 

2 
2 

14,20,22 
16,71,17 

11,18,13 
32,162,70 

9,16,13 

9,15,13 

10,13,13 

C  &  L 

4 

41,89,45 

23,61,40 

19,30,57 

19,30,57 

19,17,52 

Beale 

0 

15,47,17 

16,43,20 



Branin 
F&  R 

1 
1 

1,3,7 
11,20,13 

2,3,6 

11,21,14 

2.3,3 

2,3,3 

3,3,3 

,,,,,,, 

Bard 

1 

7,19,8 

6,17,7 

8,10,9 

8,10,9 

10,9,10 

K&  O 

0 

8,14,10 

10,13,13 

10,13,13 

10,8,16 

Osborne  1 

0 

7,33,9 

14,49,23 

49,58,229 

16,20,21 

20,19,24 

Engvall  1 

4 

2,128,40 

17,43,125 

17,41,91 

Zangwill  2 

2 

2,5,6 

3,5,4 

2,5,3 

2,5,3 

6,5,4 

Madsen 

1 

6,9,7 

6,9,10 

6,9,10 

6,9,10 

8,9,9 
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Table  5.4 


Minimum  /^c  problems 

-  Version  2 

Problem 

nz 

NPSOL 

IQP 

MQP 

PQP 

EQP 

Woods 

0 

32,88,83 

5,13,6 

7,13,7 

7,11,7 

9,10,9 

Singular 

1 

31,51,33 

31,15,34 

12,37,13 

14,37,49 

14,29,48 

Rosenbrock 

1 

24,29,38 

7,12,10 

7,11,11 

7,8,11 

7,7,11 

Chebyquad 

5 

28,57,33 

6,14,8 

6,14,8 

6,10,8 

7,10,15 

Helix 

2 

14,20,22 

9,22,13 

9,21,13 

9,14,13 

9,12,13 

Large  Res 

2 

16,71,17 

42,311,152 

C&L 

4 

41,89,45 

21,41,58 

18,31,53 

17,24,53 

16,17,50 

Beale 

0 

15,47,17 

11,29,17 



Branin 

1 

1,3,7 

2,5,4 

2,3,3 

2,2,3 

2,2,3 

F&  R 

1 

11,20,13 

11,24,17 

Bard 

1 

7,19,8 

6,28,7 

8,13,9 

8,9,9 

9,8,10 

K&  O 

0 

8,16,10 

8,12,12 

7,10,11 

9,7,16 

Osborne  1 

0 

7,33,9 

20,60,39 

16,23,21 

16,19,21 

19,18,24 

Engvall  1 

4 

14,46,24 

15,43,91 

Zangwill  2 

2 

2,5,6 

2,5,3 

2,5,3 

2,4,3 

4,4,3 

Madsen 

1 

6,9,7 

7,12,9 

7,12,9 

6,8,7 

8,9,10 
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Table  5.5 


Min 

imum  /)  problems 

Problem 

nz 

NPSOL 

IQP 

MQP 

PQP 

EQP 

Woods 

1 

6,20,14 

3,16,4 

5.16,3 

5,16,3 

13,14,4 

Singula: 

0 

11,19,12 

11,12,12 

11,10.12 

11,10,12 

16,9,12 

Rosenbrock 

0 

4,8,12 

4,9.5 

4,6,5 

4,6,5 

6,5,5 

Chebyquad 

1 

12,30,16 

9,29,13 

8,25,9 

8,25,9 

22,23,25 

Helix 

0 

8,14,11 

7,10,8 

7,10,8 

7,10,8 

9,7,8 

Large  Res 

4 

37,464,48 



C  &  L 

1 

46,77,53 

19,53,20 

19,23,20 

19,23,20 

5,22,22 

Beak 

0 

11,33,13 

10,21,11 

7,12,8 

7,12,8 

8,6,8 

Branin 

0 

7,11,8 

4,5,5 

4,5,5 

4,5,5 

7,5,5 

F&  R 

2 

9,14,10 

7,12,8 

9,13,10 

9,13,10 

10,12,10 

Bard 

0 

12,92,13 

8,47,9 

15,35,14 

16,32,17 

32,33,20 

K  &  O 

0 

10,34,14 

12,52,21 

9,23,12 

11,23,15 

34,33,30 

Osborne  1 
Engvall  1 

0 

45,182,128 



8,93,13 

'•••••• 

Zangwill  2 

0 

1,7,4 

2.7,3 

2,7,3 

2,7,3 

6,7,3 

Madsen 

0 

10,7,11 

11,7,12 

11,7,12 

11,7,12 

14,7,12 
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Table  5.6  compares  the  results  of  our  methods  and  NPSOL  on  a  group  of 
general  problems.  As  above,  the  Byrd  and  Shuitz  dropping  rule  is  not  enforced  and 
a  tolerance  of  l.e  —  i  is  used  for  determining  when  a  constraint  may  be  deleted.  The 
NPSOL  merit  function  was  used  first.  For  the  problems  on  which  our  methods 
failed  because  of  using  the  NPSOL  merit  function,  two  sets  of  results  are  listed. 
The  first  corresponds  to  using  the  NPSOL  merit  function,  and  the  second  refers  to 
the  results  obtained  from  using  the  exact  penalty  merit  function  with  a  suitable 
choice  of  p  (for  HS64,  p  =  3000;  for  Manne,  p  =  10). 

Problems  5.19-5.24.  These  problems  can  be  found  in  Hock  and  Schittkowski 
(1981).  For  each  problem,  we  list  the  number  of  variables,  n,  the  number  of  bound 
constraints,  mb,  the  number  of  linear  constraints,  ml,  the  number  of  nonlinear 
constraints,  mn,  the  starting  point,  .v°,  the  solution,  x* ,  and  the  function  value  at  the 
solution,  /  . 

Problem  5.19.   HS64. 

n  =  3,     mb  =  3,     ml  =  0,     mn  =  1. 
x°  =  (1.0,  1.0,  1.0)^ 
x'  =  (108.735,  85.1261,  204. 325)^" 
/'  =  6299.84 

Problem  5.20.    HS78. 

«  =  5,     mb  =  0,     ml  =  0,     mn  =  3. 
x^  =  (-2.0,  1.5,  2.0,  -1.0,  -1.0)^ 

x'  =  (-1.71714,  1.59571,  1.82725,  -0.76364,  -0.76364)'" 
/*  =  -2.91970 
Problem  5.21.    HS84. 

Az  =  5,     mb  =  5,     ml  =  0,     mn  =  6. 

x^  =  (2.52,  2.0,  37.5,  9.25,  6.8)^ 

x'  =  (4.53743,  2.40000,  60.0000,  9.30000,  7.00000)^ 

/*  =  -5329025 
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Problem5.22.    HS104. 

n  =  8.     mb  =  8,     ml  -~  0,     mn  -  6. 

.r°  =  (6.0,  3.0,  0.4,  0.2,  6.0,  6.0,1.0,  0.5)'" 

x'  =  (5.28632,  5.25866,  0.403204,  0.981998,  9.26043. 

6.76707,  1.95311,  1.99581)'" 
/•  =  1.00000 
(Note    that    this    solution    is    different    from    the    one    given    in    Hock    and 
Shittkowski.) 

Problem  5.23.    HSlli. 

n  =  10,     mb  =   10,     ml  =  0.     mn  =  3. 

x^  =  (-2.3.    -2.3,  -2.3,  -2.3.    -2.3.  -2.3,  -2.3,  -2.3,  -2.3.  -2.3)^ 

x'  =  (-3.20231,  -1.91237,  -0.244427,  -6.56118,  -0.723098. 

-7.27423,  -3.59724,  -4.02032,  -3.28838,  -2.33437)'' 
/*  =  -47.76109 
(The  solution  listed  in  Hock  and  Shittkowski  for  this  problem  is  not  accurate. 
Our  solution  agrees  with  the  one  given  in  Nocedal  and  Overton  (1985)). 
Problem  5.24.    HSil8. 

n  =  15,     mb  =  15,     ml  =   17,     mn  =  0. 

x^  =  (20.0,  55.0,  15.0,  20.0,  60.0,  20.0,  20.0, 

60  0,  20.0,  20.0,  60.0,  20.0,  20.0,  60.0,  20.0)'" 
x'  =  (8.0.  49.0,  3.0,  1.0,  56.0,  0.0,  1.0. 

63.0,  6.0,  3.0,  70,0.  12.0.  5.0.  77.0,  18.0)'" 
/*  =  664.8204 

Problem  5.25.    Manne,    This  is  a  smaller  version  of  a  problem  presented  in 
Murtagh  and  Saunders  (1982), 

n  =  30,     mb  =  30,     ml  =  10,     mn  =10. 
x^  =  (3,05,  3,1.  3,2.  3.3,  3.4,  3.5,  3.6,  3.7,  3.8,  3.9, 
0.95.  1.0,  1.0,  1,0,  1,0.  1,0.  1,0.  1,0,  1,0,  1,0. 
0.1,  0.1,  0.1,  0.1,  0.1,  0,1,0,1.  0.1,0.1,  0. 1)'' 
x'  =  (3.05,  3.12665.  3.21445.  3.30399,  3.39523. 

3.48790,  3,58172,  3,67638,  3,77154.  3.86667, 
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0.95,  0.968400,  0.99783,  1.02818,  1.05966, 
1.09230,  1.12611,  1.16116,  1.19758,  1.21394, 
0.766504^-01,  0.877975e-01,  0.895466e-01, 
0.91237^-01,0.92616^-01, 
0.9381446^-01,  0.9466694^-01,  0.951548^-01, 
0.951273^-01,  0.116000)^ 
-2.670099 

Table  5.6 


General  Problems 

Problem 

nz 

NPSOL 

IQP 

MQP 

PQP 

EQP 

HS  64-1 

2 

34,56,51 

17,38,33 

HS64-2* 

2 

39,88,63 

41,45,227 

41,45,227 

47,49,238 

HS78 

2 

8,9,10 

6,7,8 

7,8,9 

7,8,9 

7,8,9 

HS84 

0 

5,9,6 

2,4,3 

2,4,3 

2,4,3 

5,4,6 

HS104 

7 

7,11,8 

6,6,7 

5,6,6, 

5,6,6 

6,8,11 

HSlll 

7 

56,57,59 

56,57,60 

56,57,60 

56,57,60 

56,57,60 

HS118 

0 

14,33,30 

9,29,11 

9,25,11 

9,23,11 

17,17,17 

Manne-1 

7 

12,37,13 

9,31,12 

Manne-2* 

7 

9,31,12 

40,64,41 

40,64,41 

58,67,110 

*  using  the  exact  penalty  merit  function 

The  results  demonstrate  that  the  NPSOL  merit  function  can  generally  be  used 
for  all  our  methods  when  the  IQP  method  performs  only  one  minor  iteration  per 
major  iteration.  In  that  case,  our  other  methods  perform  very  similarly  to  the  IQP 
method,  and  the  QP  multiplier  estimates  should  have  the  correct  signs.  On  the 
other  hand,  for  problems  on  which  the  IQP  method  takes  many  minor  iterations 
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during  the  course  of  one  major  iteration,  the  other  methods,  which  are  forced  to 
terminate  the  subproblem  prematurel)',  suffer  from  inaccurate  multiplier  estimates, 
(n  most  cases,  the  exact  penalty  merit  function  can  be  used,  although,  as  mentioned 
above,  it  is  extremely  sensitive  to  the  proper  choice  of  p.  For  HS64,  setting 
p  =  3000  provided  satisfactory  results;  however  for  p  =  2280,  the  algorithms  failed. 

The  problems  for  which  nz  is  relatively  large,  HS104,  HSlll  and  Marine,  were 
chosen  in  order  to  compare  our  methods  with  NPSOL  on  problems  for  which  the 
dimension  of  the  projected  Hessian  matrix  is  not  significantly  smaller  than  the 
diracnsion  of  the  full  Hessian.  Problems  HS78  and  HSlll  are  equahty  constrained 
problems.  Therefore,  the  active  set  is  immediately  identified  and  one  would  expect 
the  projected  methods  to  perform  well,  even  when  nz  is  large.  Problem  Manne  is  a 
problem  for  which  constraints  must  be  added  to  the  working  set,  but  are  later 
deleted  from  it  Our  non-IQP  methods  suffer  from  this  need  to  often  revise  the 
working  set.  Problem  HS118  is  a  linearly-constrained  problem,  for  which  nz  is  zero 
at  the  solution,  but  is  14  at  the  first  iteration.  This  indicates  that  many  constraints 
must  be  added  to  the  working  set.  Since  the  EQP  method  can  add  only  one 
constraint  per  iteration,  it  is  not  surprising  that  it  takes  ainiost  twice  as  many  major 
iterations  as  our  other  methods.  The  other  methods,  however,  take  many  more 
minor  iterations. 

Problem  5.26.  Penalty  Function  I.  This  problem  can  be  found  in  More, 
Garbow  and  Hillstrom  (1981). 

n  =  50,     mb  =  50,    ml  =  0,     mn  =  0. 


0  _ 


(1,  -1. 1,1,     ■  •  .  1,  -n'^ 


x'  =  (0.71621e-01,  0.71621e-01,    •  •  •  ,  0.71621^-01)^ 
/*  =  0.4313635^-01 
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Table  5.7 


Penalty  Function  I 

delete  tol 

Method 

9.99e+5 

9.99e-2 

NPSOL 

9,131,27 

9,131,27 

IQP 

16,164,44 

16,164,44 

MQP 

9,59,27 

>200 

PQP 

9,59,27 

>200 

EQP 

9,10,27 

29,30,93 

This  problem  is  an  example   where  an  EQP  method  works  best.    On  this 
problem  the  IQP  method  picks  up  all  50  bound  constraints  at  the  first  iteration  and 
then  must  subsequently  delete  them.    Smce  the  EQP  method  cannot  pick  up  many 
constraints  at  once,  it  does  not  go  through  this  process  of  unnecessarily  modifying 
the  working  set.    Nevertheless,  the  IQP  method  does  recover  from  its  initially  bad 
start.   The  MQP  and  PQP  methods,  which  also  add  many  constraints  during  the  first 
Iteration,   can   recover   faster   than   the   IQP   method   provided   that   a   very    loose 
dropping  tolerance  is  used.    If  any  reasonably  small  tolerance  level  is  used,  these 
methods  suffer  from  not  being  able  to  delete  constraints  until  much  effort  has  been 
spent  in  finding  the  minimum  on  the  wrong  subspace.    We  note  that  the  MQP  and 
PQP  methods  produce  identical  iterates  for  this  problem.   This  is  due  to  the  fact  that 
at  the  first  iteration  the  IQP  adds  many  constraints  to  the  working  set.  but  does  not 
delete  any  of  them  until  subsequent  iterations.    Therefore,  at  the  end  of  the  first 
subproblem  all  of  the  non-EQP  methods  are  at  the  same  point.   For  subsequent  QP 
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subproblems,  both  the  MQP  and  the  PQP  methods  start  with  the  same  initial 
working  sets.  Since  constraints  are  no  longer  being  added  within  the  subproblems, 
the  MQP  method  cannot  drop  any  constraints  within  the  QP,  and  thus  is  the  same  as 
the  PQP  method  for  this  problem. 


5.2  Conclusions 

The  results  indicate  that  our  IQP  method  is  the  most  robust  of  our  methods.  In 
many  cases  the  superiority  of  the  IQP  method  is  due  to  the  ability  to  use  the  NPSOL 
merit  function  rather  than  the  exact  penalty  merit  function.  However,  even  for  the 
minimax  and  minimum  /j  problems,  for  which  a  merit  function  is  available  for  all 
the  methods,  the  IQP  method  performs  the  best.  Yet,  for  problems  on  which  none 
of  the  methods  have  any  difficulty,  the  IQP  method  may  be  slightly  worse  due  to  the 
fact  that  it  often  takes  more  minor  iterations. 

The  MQP  and  PQP  methods  almost  always  perform  identically.  These 
methods  can  differ  only  in  a  case  where  the  IQP  would  pick  up  many  constraints 
and  then  drop  some.  Usually,  if  there  is  a  difference  between  the  two  methods,  it  is 
evidenced  at  the  first  iteration.  During  subsequent  iterations,  the  IQP  may  modify 
the  working  set  by  adding  new  constraints,  or  deleting  old  ones,  but  it  rarely  adds 
many  constraints  and  then  drops  some  within  one  QP  once  the  computation  has 
passed  the  first  few  iterations.  Therefore,  at  that  stage  the  MQP  and  PQP  methods 
are  identical. 

Since  the  MQP  and  PQP  methods  can  fail  when  used  with  the  NPSOL  merit 
function,  a  possible  improvement  would  be  to  use  a  merit  function  which  is  more 
suitable  to  these  methods.  It  would  appear  that  such  a  merit  function  might  be 
defined  as  an  adaptive  method  which  does  not  reh  on  multiplier  estimates  when  far 
from  the  solution,  but  which  allows  for  fast  local  convergence  by  considering  the 
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multipliers  once  the  working  set  has  settled  down. 

The  EQP  method  is  preferable  in  cases  in  which  the  IQP  picks  up  many 
constraints  but  then  must  delete  them  all  from  the  working  set.  Since  the  EQP  can 
pick  up  only  one  constraint  at  a  time,  it  does  not  have  the  chance  to  add  many 
constraints  that  are  incorrectly  determined  to  be  active.  However,  for  a  problem  in 
which  many  constraints  must  be  made  active,  the  EQP  will  take  many  more  steps 
(and  thus  will  make  more  function  calls)  in  order  to  add  each  constraint  separately 
to  the  working  set. 

In  general,  our  results  indicate  that  the  projected  Hessian  updating  methods  are 
most  successful  when  the  dimension  of  the  projected  Hessian  is  small  and  the  active 
constraints  are  correctly  identified  in  the  early  stages  of  the  computation.  These 
methods  perform  well  even  when  the  dimension  of  the  projected  space  is  large, 
provided  that  the  active  set  is  easily  identified.  Performance  of  these  methods  is 
worst  on  problems  for  which  the  working  set  changes  many  times,  since  the  method 
does  not  get  a  chance  to  build  up  second  derivative  information  on  the  projected 
subspace. 

The  most  important  question  is  whether  projected  Hessian  updating  methods 
have  a  useful  role  to  play  in  future  choices  of  the  best  methods  for  solving  nonlinear 
programming  problems.  We  feel  that  our  results  for  the  IQP  algorithm  are 
encouraging  in  that  regard.  Although  the  NPSOL  results  are  very  good,  our  IQP 
algorithm,  using  the  NPSOL  merit  function,  can  give  somewhat  better  results,  as 
demonstrated  by  Table  5.6.  Note  that  essentially  the  only  difference  between  the 
two  methods  is  that  we  recur  the  projected  Hessian  rather  than  the  full  second 
derivative  matrix.  The  numerical  results  indicate  that  projected  Hessian  updating 
methods  may  be  useful  in  practice. 
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Appendix 
An  Image-Processing  Problem 


Introduction 

We  consider  the  following  problem  which  arises  in  image  processing 
applications  (Schwartz,  private  communication).  The  matrices  A,B  and  C  are 
positive  definite. 

min  f(u,v,w)  =  u^Au  +  v^fiv  +  w^Cw  (A.l) 

subject  to        u^u  =  \ 


w   w  =    1 

w^■  =  0 

K^H    =   0 
rV  =  0 
The  solution  is  not  unique,  since  any  of  the  vectors  //,  v  and  i>  can  be  replaced  by  its 

negative.    We  will  denote  the  ;th  constraint  by  c,(;<.\',uj.    Let  g  denote  the  gradient 

of  th,    objective  function  and  let  A  be  the  6x9  matrix  whose  rows  consist  of  the 

transposes  of  the  gradients  of  the  constraints.    Let  G  and  G,  denote  the  Hessians  of 

the  objective  and  constraint  functions,  respectively.    We  have  then 


2Au 
IBv 
ICw 
and 


2«    0      0     V   u'    0 
0    2v     0     «    0    w 
0     0    Iw   0    It    V 
We  now  define  Z  to  be  a  matrix  with  orthogonal  columns  spanning  the  null  space  of 
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A.   When  the  constraints  are  satisfied,  such  a  Z  is  given  by 


Z  = 


V 

0 

H' 

—  u 

w 

0 

0 

—  V 

—  u 

Since  A^X  =  5,  we  have 

u^Au 

w^Cw 
(v^A«  +  uJBv) 
{w'^Au  +  u^Cw) 
(^w'^Bv  +  v^Cw) 

and  the  projected  Hessian  of  the  Lagrangian  function  has  the  form: 


X  = 


Iv'^Av-luJAu  +  lu'^Bu-lv^Bv  w^Au+ u^Cw -lu^Bw 


-  wtBv-  v^Cw  +  Iv^Aw 


^Tau  +  u^C^-Iu^Bw  2w^flw-2v^flv  +  2v^Cv-2w^C>v  -  v^Au- u^Sv  + 2v^Cu 

-w^flv-v^Cw  +  2w^Av  -v^Au-u^flv  +  2u^Cv  2>v^A>v  -  2u^A«  +  2«^Cu- 2w  Cw 


Description  of  the  Algorithm 

The  problem  under  consideration  is  a  special  case  of  the  general  nonlinear 
programming  problem  (NP)  discussed  in  this  thesis.  Note  that  all  the  constraints 
are  equalities,  so  all  of  the  SQP  methods  discussed  in  this  thesis  reduce  to  solving  an 
EQP  at  each  step.  Since  in  this  application  the  second  derivative  matrices  {G,}  are 
readily  available,  there  is  no  advantage  to  using  quasi-Newton  updating;  instead  the 
algorithms  can  form  Z^WZ  directly. 

A  special  feature  of  Problem  (A.l)  is  that  it  is  relatively  easy  to  ensure  that  the 
constraints  remain  satisfied.  If  the  new  iterate  does  not  satisfy  the  constraints,  the 
vectors  u,v,  and  w  can  be  orthonormalized  by  applying  the  Gram-Schmidt  process. 
Thus  c  will  always  be  zero,  and  the  p,  term  defined  by  (1.10a)  can  be  ignored.  The 
Newton  step  is  thus  defined  simply  by 
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Z^WZp,  =  -Z^g 

P    =  Zp: 

We  present  an  SQP  method  for  solving  (A.l),  in  which  we  take  advantage  of 
the  pecial  structure  of  the  constraints  The  method  we  propose  generates  a 
sequtFice  of  feasible  points  which  converge  to  a  first-order  minimizer  of  (A  1).  If 
the  projected  Hessian  'S  not  positive  definite  at  the  solution,  we  introduce  a 
perturbation  in  order  to  move  away  from  a  saddle  point  However,  since  the 
problem  r/ay  not  have  a  strong  local  minimum,  this  is  done  only  a  few  times. 

START  with;/,  r'\  kv°. 
We  will  use  the  notation  (u.v.w)  to  denote 

Set  count  =  0. 

evaluate  Z^  Z*'^S  X*,  Z'^^W'^Z''. 

WHILE  (  |iZ*'^g''||>€  or  Z'^^w'^z''  is  not  positive  definite  and  coiint<4  ) 
DO 
evaluate  Z*.  Z*'?^  \*,  Z'''w''Z''. 

compute  the  modified  Cholesky  factorization  of  Z^'^W'^Z*. 
(this  in  effect  produces  a  factorization  of  Z*'lV*Z*  +  £* .  where  E* 
is  a  diagonal  matrix  which  is  added  so  that  the  resulting  matrix  is 
positive  definite  (see  Gill,  et.  al.  (  1981).) 
IF  Z'''w''Z''  is  not  positive  definite 
THEN  set  count  =  count+\ 
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IF  ||Z*^/||>€ 

THEN  set  the  r.h.s.  to  Z*'^* 

ELSE  set  the  r.h.s.  to  a  random  positive  vector 
solve  (Z*V*Z*  +  £*)/7*  =  -r.h.s. 
Setp*  =  ZV* 
Set  a  =  1 

REPEAT 

Set  (m"',v*,w'')  =  («*,v*,w*)  +  ap* 

use  Gram-Schmidt  to  orthonomalize  (m*,v''',w'*") 

Set  a  =  a/2 
UNTIL  (  /(w*,v'+,vv^)</(i/,v*,H''^)  ) 
Set  a*  =  a 

Set(t/**\v*^i,vv*^i)  =  (,/*,v^,w*) 
Set  A:  =/t+l 
evaluate  Z*.  Z*'g*,  \*,  Z'^V^^Z*. 

ENDWHILE 

In  the  context  of  image  processing,  Problem  (A.l)   arises  naturally  in  three 
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dimensions.  Note,  however,  that  the  algorithm  can  be  used  to  solve  the  general 
problerr;  in  an\  number  of  dimensions     The  general  problem  is 

min    f(xi,X2.  ■      ■  ,x„)  =  x/a^x^  +  xi^Aix-  +    ■  ■  ■    +  xJa„x„ 

St.    .r/,c,  =   1 

Xi^Xj  =  0        i^; 
and  the  matrices  .4,  are  all  positive  definite.    The  problem  then  has  n'^  unknowns  and 

^^ constraints.    The  matri.x  Z.  whose  columns  form  a  basis  for  the  null  space 

2 

2 

of  the   constraint   gradients,  has    — - —  columns.  The   algorithm  therefore  solves 

2  2_ 

X  — - —  systems  of  equations  at  each  iteration  A  nice  feature  of  the  two- 
dimensional  case  IS  that  the  iNewton  step  satisfies  the  orthogonality  constraint,  so 
that  the  new  iterate  can  be  found  simply  by  normalizing  the  Newton  step  In  the 
following,  we  will  restrict  our  attention  to  the  original  problem  as  posed  in  R". 
Clearly,  a  similar  analysis  can  be  used  to  prove  convergence  in  R". 


Global  Convergence 

The  algorithm  generates  a  sequence  of  feasible  points  Thus  to  prove 
convergence  to  a  first-order  minimizer  of  Problem  (1).  it  is  sufficient  to  show  that 
linillZ'^'"^*!!  -  0.    We  show  that  we  can  choose  a  direction.  /7*,  from  the  current  point 

(m*,\'',w*)  along  which  the  objective  function  decreases  initially.  We  travel  along 
thi^  lirection  to  a  new  point  {u~  .v^  .w~)  and  apply  Gram-Schmidt  to  (u'^  ,v~  ,\v*)  to 
genjfate  (w**^v**\vv*^-)  which  is  a  feasible  point.  The  step  along  the  direction  p* 
can  be  chosen  so  that /(w*"',v*'^u*"')  < /(;/,v•^i^•S.  We  show  that  the  step  to 
the  constraints  does  not  impede  descent  and  that  the  step  length,  a*^  does  not 
become  arbitrarily  small.    Finally,  we  show  that  in  the  limit,  the  projected  gradient 
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tends to  zero  and  thus  the  algorithm  is  globally  convergent. 

Any  direction  along  which  the  objective  funtion  decreases  initially  is  called  a 
descent  direction.  If  we  consider  the  local  Taylor  series  approximation  to  /,  it  is 
clear  that  if  ^*^/7*<0,  then  /?*  is  a  descent  direction. 

Lemma  1:  p*  is  a  descent  direction  for  /. 

Proof:  By  definition,  (Z*V*Z*+ £*);?*  =  -Z*'f*.   Thus 

Since  Z*'^H'*Z*+£*  is  positive  definite,  ^*p*^  must  be  negative  if  p*?^0,  i.e.,  p'^=^0.  □ 
Let  p*  =  (a,b,c)^  and  p*  =  max(|a|,|fe|,|c|). 
Lemma  2:  At  (u'*' ,v*  ,w*)  the  constraint  violations  are  0(a*V*V 
Proof: 

v^^*  =  1  +  a''a^  +  a'^'b^ 
w"V^  =  1  +  a*'i,2  +  ^k^^l 

w     V     =  -a  (3  4-  a  a  +  a^c 

V*  w"^  =  -OL'b  +  a*i>  +  a*^<3c 

w*  w"^  =  a*c  -  a*^"bc  +  a*c 
Lemma  3:  The  step  from  (i<'^,v*,iv*)  to  the  constraints,  i.e.,  the  change  given 
by  the  Gram-Schmidt  process,  is  0(a*^p*^). 
Proof:  The  new  iterates  are  given  by 

\\u    \\ 

||v"    -    («**l'v^)M**l|| 


(M^^^)«^ 


V      ~  z 
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.*+ir,.,*-..>^l  /,,*+ir,.,  +  x,,«r*l 


The  size  uf  the  step  taken  to  the  constraints  is  thus 

I  11"  H  II"  II  ) 


=     I  ;j        ;y       —    7  ly        


!<        J/       —    J.U 


-    1 


=  Vn: 


(l-||«''|i)^  =  6>(a*"p'-")    by  Lemma  2. 


ii  =  ^v-^v 


-  -  2v^V^i  +  1 


,-         {u*\-)u' 

_r  _ 

-  2v-^ 

.r    _ 

.r    .         2(u'  v~)           (u*  v~)  u*  u~ 

V      V ^ <—  ->■   -^ '■ — ^ 

\                      u'  u'                    u~  u' 

_r  ^ 

u     u 

.r 

»        U                                  U         .' 

■,  *  -^ 

-    -r    _  ,_M       V 


+r    +  _    .r    _  _ 

V     v     —  lev     V     —  2c 


where 


(  T      +'■     +  /      +'     '^\2 

^T  +         2u     V             (u     V    )^ 
"■     V     —  : +  


Thus,  c  =  0(  1  +■  ya*'p*')  .    It  follows  that 


-  v^ll  =    V(i  _  ||,.  +  ||)2  ^  0(a'^p'^)  =  0(ot*V*')  as  abov< 
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Similarly,  it  can  be  shown  that 


;*"!  -  w^ll  =  V,v  +  7-v^+  -  2vv^'"w*^l  +   1 


V^ 


||.v*||)2+0(aV)  =  0«.V) 


Given  that  the  step  to  the  constraints  is  0(a*  p*"),  it  follows  that 

/(w*^\v*^\h'*+1)  =/(  (M  +  ,v^.w  +  )  +  O(aV')  ) 

=  /(  (M*,v*,vv*)  +  aV*  +  «*V*'  ) 

=  /(M*,v*,w*)  +  a*^*p*  +  0(aV'). 
The  method  of  computing  the  factorization  guarantees  that  Z*  W*Z*  +  £*  is  bounded 

away  from  singularity.    Let  t]*  and  5*  be  the  minimum  and  maximum  eigenvalues  of 

2*^*^2*^  +  £*,  respectively.    It  follows  that  ti''||Z*'^*||  ^  |[p*||  ^  C*||Z*'^*||.  Thus, 

from  the  proof  of  Lemma  1  it  follows  that  a*^*p*  is  a  negative  term  of  order  a*p*^ 

Therefore,  there  must  exist  an  a*  such  that /(w*'^\v''^^\h'*"^^)</(w'^,v*,w*). 

Lemma  4:  a*  is  bounded  away  from  zero. 

Proof:  Since  the  method  generates  only  feasible  points,  ||i<*||=  ||v*||=  ||vv*||=l. 
As  noted  above,  \\pr\\  is  bounded.  Moreover,  we  note  that  since 
^*  =  (2/l;/,  2S^■^  2Ch'*)^,    ||^*||  <  2V3A/.    Thus,    p*    is    bounded.     Clearly,    if   we 

choose  a    <  — ,  where  K  is  an  appropriate  constant,  we  obtain  a  decrease  in  /. 

Since  p*  cannot  become  arbitrarily  large,  a*  is  bounded  away  from  zero.    Note  that 
when  p*  is  sufficiently  small,  we  can  take  a*=  1. 

Theorem:  If  (u  ,v  ,h'  )  satisfy  the  constraints,  then  either  the  sequence 
(u  ,v  ,w  )  generated  by  the  method  terminates  at  a  stationary  point  of  Problem  (1) 


liminf  ||Z*'^*||  =  0. 

k~cc 

Proof:  Assume  that  the  conclusion  of  the  theorem  is  false.   Then  there  exists  an 
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6>0  buch  that  ||Z*'^^*||  ^  €  for  all  k.    We  have  that 


/(«*,v*V)  -/(//  +  \v'^^\i/"^)  =  -(aVp     +  0(aV*))- 


Since  /  is  bounded  below,  by  taking  the  sums  of  both  sides  of  the  equation,  we 
obtain,  that 

k=0  k=0 

Since  the  assumption  is  that  ||Z*g*||  ^  e  >  0  for  all  k.  and  ■ti*>0  is  guaranteed  by 
the  method  of  computing  the  modified  Cholesky  decomposition,  we  have 

lim  a*  =  0. 
However,  this  contradicts  Lemma  4.    Hence  the  assumption  was  false. 
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